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ABSTRACT 


The  program  SAGA  exploits  short  arc  orbital  constraints  i n  effecting  the 
adjustment  of  observations  made  by  geodetic  tracking  nets  embracing  both  optical  systems 
(e.g.,  PC-1000,  MOTS)  and  electronic  ranging  systems  (Lasers,  Secot ,  Geoceiver,  etc.) 
Provisions  are  made  for  consideration  of: 

a)  random  errors  in  the  observations  and  in  the  iiming  of  observations; 

b)  serially  correlated  errors  in  observations; 

c)  errors  in  the  adopted  location  of  the  center  of  mass; 

d)  systematic  errors  governed  by  erior  mouels  having  coefficients 
subject  to  a  priori  contrail, 

The  overall  tracking  net  can  include  an  indefinitely  large  number  of  stations  (many 
hundreds)  o*,  long  as  no  more  than  fifteen  participate  successfully  in  the  observations  of 
any  pass.  All  orbital  state  vectors  are  treated  as  unknown  and  no  limits  are  sot  on  the 
number  of  state  vectors  thct  can  be  solved  for  simultaneously.  Allowances  are  made  in 
optical  error  modelling  for  reinitialization  of  error  coefficients  that  becomes  necessary 
when  any  station  exposes  more  than  one  plate  on  a  given  pass.  In  the  case  of  electronic 
tracking,  up  to  three  dropouts  in  tracking  can  be  accommodated  for  each  station  on 
each  pass  with  appropriate  reinitialization  of  error  coefficients.  A  maximum  of  over 
250  error  coefficients  can  be  exercised  in  the  reduction  of  each  pass.  This  becomes 
computationally  feasible  by  virtue  of  algorithms  providing  the  solution  to  problems  in 
what  is  termed  second  order  partitioned  regression. 
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INTRODUCTION 


In  previous  series  of  investigations  conducted  by  DBA  Systems  for  AFCRL 
(Brown,  Bush,  Sibol  (1963),  (1964);  Brown  (1964b),  we  developed  most  of  the  theoretical 
framework  for  the  present  undertaking,  namely  the  development  of  a  general  and  advanced 
computer  program  for  short  arc  geodesy.  This  program,  hereafter  referred  to  as  SAGA  (Short 
Arc  Geodetic  Adjustment),  has  drawn  on  and  benefited  from  experience  gained  from  three 
preceding  programs  that  were  also  based  on  the  theoretical  development  referred  to  above. 
These  predecessors  consist  of  the  following  programs  developed  by  DBA: 

(1)  GDAP  (GEOS  Data  Adjustment  Program)  developed  for  NASA  Goddard 
during  the  period  1966-1967  for  the  short  arc  reduction  of  observations 
of  geodetic  satellites  (Lynn,  1967); 

(2)  MCT  (Method  of  Continuous  Traces)  developed  for  AFCRL  during  the 
period  1966-1967  for  the  recovery  of  geodetic  positions  from  measure¬ 
ments  of  sun-illuminated  passive  satellites  recorded  against  the  stellar 
background  (Brown,  Trotter,  1967); 

(3)  NAP  (Tracking  Network  Analysis  Program)  developed  for  NASA  Goddard 
during  the  period  1967-196?  for  long  arc  orbital  reduction,  terrestrial, 
lunar,  or  interplanetary,  with  options  for  recovery  of  station  coordinates 
end  error  model  coefficients  (Lynn,  et  al,  1969). 

Properties  of  SAGA  relative  to  GDAP,  MCT  and  NAP  are  indicated  in  Table  1.  All  four 
programs  are  capable  of  exercising  short  arc  orbital  constraints  in  an  unlimited  multi-epoch 
mode  (that  is,  the  number  of  state  vectors  that  can  be  solved  for  simultaneously  is  without 
present  limit) .  In  addition,  NAP  can  exercise  long  arc  constraints  and  can  accommodate 
extraterrestrial  orbits  (it  employs  a  general  n  body  integrator  capable  of  integrating  through 
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Short  arc  constraints  (unlimited  no.  epochs) 
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Long  arc  constraints  (unlimited  no.  epochs) 
Extraterrestrial  orbits 


Geometric  option  (orbital  constraints  not  used) 
Recovery  of  station  locations 

(a)  limited  number 

(b)  unlimited  number 

Recovery  of  coefficients  of  potential  function 
Recovery  of  coordinates  of  center  of  mass 
Recovery  of  tracking  error  coefficients 

(a)  optionally  reinitialized  after  every  pass 

(b)  optionally  stable  over  specified  sets  of  passes 

Consideration  of  random  timing  ciror 

Option  for  consideration  or  :eiiully  correlated  errors 
by  Autoregressive  Fei.-di.ack 

Solution  of  genera!  normal  equations  by  First  Order 
Partitioned  Regression 

Solution  of  general  normal  equations  by  Second 
Order  Partitioned  Regression 

Tracking  systems  accommodated 
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(a)  Optical  (PC-1000,  MOTS,  Baker-Nunn, 
BC-4,  etc.) 

(b)  Electronic  Ranging  (Lasers,  SECOR, 

GRARR,  Radars) 

(c)  Microwave  Interferometer  (MINITRACK) 

(<0  Noncumulative,  one  way  doppler  (TRANET) 
(e)  Cumulative,  one  way  doppler  (GEOCEIVER) 


X 

X 

X 

X 


Structure  of  program  and  organization  of  inpuf/output 
optimized  primarily  for 

(o)  Recovery  of  station  locations 
(b)  Recovery  of  precise  orbits 
.  (c)  Recovery  of  tracking  error  coefficient* 
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15.  Program  operational  on 


(a)  CDC  3100/3200  (16K  core,  4  mag.  tape  units) 

(b)  CDC  3300  (Naval  Research  Lab.) 

(c)  UNIVAC  1230  (NASA  Goddard) 

(d)  IBM  360/75/91/95  (NASA  Goddard) 

(e)  IBM  7044/7094  (AFCRL) 
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Extension  underway  to  add  indicated  capability. 
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spheres  of  influence  while  allowing  the  gravitational  potential  of  the  currently  dominant 
center  of  attraction  to  be  represented  by  solid  spherical  harmonics  up  to  degree  and  order 
(n,m)  =  (24,24).  While  NAP  can  be  employed  for  satellite  geode---  its  primary  intent  and 
organization  are  directed  toward  precise  determination  of  orbits  with  appropriate  consideration 
being  given  to  (a)  serially  correlated  errors  of  fracking  systems,  (b)  systematic  errors  of 
tracking  systems,  (c)  errors  in  locations  of  tracking  stations  and  (d)  errors  In  coefficients  of 
the  potential  function  (at  this  writing  the  capability  (d)  is  In  the  process  of  being  implemented 
in  NAP).  Two  advanced  features  first  introduced  in  NAP  have  been  incorporated  Into  SAGA, 
namely,  solution  of  the  normal  equations  by  means  of  second  order  partitioned  regression 
and  consideration  of  serially  correlated  errors  by  auto  rag  resslve  feedback.  In  Sections  4 
and  5  of  this  report  we  shall  provide  the  detailed  development  of  the  theories  of  partitioned 
regression  and  autoregressive  feedback. 

All  four  programs  can  recover  coordinates  of  tracking  stations.  However,  in 
forming  the  reduced  system  of  normal  equations ,GDAP  and  NAP  retain  the  system  in  core. 

This  places  a  definite  limit  to  the  number  of  stations  that  con  be  solved  for  simultaneously 
(typically  to  20  to  60  depending  on  computer).  On  the  other  hand,  SAGA  employs  the 
logical  development  originally  proved  in  MCT  wherein  the  reduced  normal  equations  are 
generated  piecewise  in  core  but  are  cumulatively  formed  on  an  external  file  (magnetic  tape 
or  disk) .  By  this  means  it  becomes  practical  to  accommodate  an  overall  tracking  network 
embracing  literally  hundreds  of  unknown  stations.  The  primary  restriction  is  that  only  a 
limited  subset  of  stations  Is  regarded  as  participating  successfully  in  the  observation  of  any 
given  pass  (in  SAGA  the  number  is  limited  to  a  maximum  of  fifteen).  Such  a  restriction  is 
of  no  practical  consequence  in  actual  short  arc  operations,  for  rarely  would  as  many  as 
fifteen  stations  participate  on  a  given  pass,  mud-  less  all  be  successful. 

Error  model  coefficients  appropriate  to  each  channel  of  observations  can  be 
carried  as  adjustable  constrained  parameters  in  all  four  programs.  In  MCT  and  SAGA,  all 
exercised  error  coefficients  are  regarded  as  unstable  from  pass  to  pass  and  thus  are  auto¬ 
matically  reinitialized  on  each  pass.  GDAP  and  NAP  also  have  this  capability  but  art 
somewhat  more  flexible  in  that  any  desired  subset  of  error  coefficients  can,  on  option,  be 
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treated  as  stable  over  a  specified  set  of  posses.  For  example,  timing  bias  from  a  given 
station  can  be  regarded,  when  desired,  as  stable  over  certain  specified  passes,  rather 
than  being  reinitialized  on  each  pass  as  are  the  other  coefficients.  This  capability  was 
not  incorporated  into  SAGA  primarily  because  it  complicates  considerably  the  logic  and 
set  up  of  the  program  and  has  proven  to  bo  a  feature  that  is  not  often  exercised  in  practice. 

Of  the  four  programs  only  GDAP  has  the  option  to  perform  geodetic  reductions 
in  a  strictly  geometric  mode.  This  option  was  not  incorporated  Into  the  other  programs 
because  experience  with  GDAP  demonstrated  the  clear  superiority  in  satellite  geodesy 
of  the  short  arc  mode  over  the  geometrical  mode,  in  particular,  recovery  of  tracking  error 
coefficients  has  been  found  to  be  far  superior  in  the  short  arc  mod**.  Comparative  analyses 
of  the  short  arc  versus  the  geometrical  op;  -each  are  made  in  Brown  (1967a)  and  Brown  (1968). 

SAGA  incorporate  cvs.ain  unique  capabilities  which  experience  with  the  other 
programs  indicated  would  be  desirable.  One  is  consideration  of  random  timing  error.  This 
was  included  primarily  to  make  proper  allowances  in  the  reduction  of  PC-1000  chopping 
observations  of  passive  satellites  in  view  of  studies  indicating  an  rms  mechanical  jitter  of 
about  0.3  ms  in  the  operation  of  the  chopping  shutter.  Inasmuch  as  relatively  close 
satellite  passes  can  cross  the  plate  of  a  PC-1000  at  rates  up  to  lQmm/sec.,  an  rms  error  in 
timing  of  0.3  ms  can  be  equivalent  to  as  much  as  3  micr  - •>>  on  the  plate  and  thus  be 
comparable  to  plate  measuring  errors.  Although  it  remains  yet  to  be  determined,  random 
timing  error  might  also  potentially  be  of  significance  with  the  Geocelver,  particularly  if 
rms  noise  levels  in  phase  determination  amounting  to  only  about  0.1  m  are  achieved  as 
projected  in  the  design.  Inasmuch  as  range  rate  of  an  observed  satellite  can  amount  to  as 
much  as  5000rrvAec,  an  rms  error  in  timing  of  as  little  as  20  microseconds  can  be  equivalent 
to  the  expected  rms  error  of  0.1  m  in  phase  determination.  Therefore,  to  be  on  the  safe  side, 
SAGA  also  makes  provisions  for  the  possibility  of  significant  random  timing  error  in  electronic 
observations. 

Another  uniq'  e  feature  of  SAGA  is  its  ability  to  take  rigorously  into  account 
errors  in  the  adopted  center  of  mass  of  the  earth.  As  was  pointed  out  in  Brown  (1967),  when 
one  elects  in  a  short  arc  reduction  to  hold  fixed  the  coordinates  of  a  selected  station,  one 


thereby  implicitly  defines  what  is  to  be  regarded  os  the  location  of  the  Earth's  center 
of  mass.  The  error  in  this  implicitly  adopted  center  of  moss  does  have  an  effect  (albeit 
a  rather  weak  effect  in  most  cases)  on  recovered  coordinates  of  tracking  stations.  In  the 
case  of  optical  tracking  networks  of  continental  extent,  the  effect  of  such  errors  is  very 
nearly  equivalent  to  that  of  an  error  in  scale*  Through  a  series  of  exercises  conducted 
with  GDAP,  Lynn,  (1969)  established,  for  example,  that  an  error  of  50  meters  in 
the  vertical  component  of  the  adopted  location  of  the  center  of  mass  (verticc  i,  that  It, 
with  respect  to  the  fixed  station)  can  cause  an  error  In  scale  of  about  1:200,000  in  a 
continental  network.  In  view  of  such  findings,  we  incorporated  into  SAGA  the  capability 
of  treating  the  coordinates  of  the  adopted  center  of  mass  as  constrained  parameters.  Tbit 
means  that  in  sufficiently  strong  tracking  networks  of  continental  extent,  the  possibility 
emerges  of  improving  the  location  of  the  center  of  mass  relative  to  the  origin  of  the 
adopted  datum.  In  weaker,  more  limited  networks  the  main  benefit  of  carrying  coordinates 
of  the  center  of  mass  as  constrained  parameters  lies  in  the  more  comprehensive  and  realistic 
error  propagation  that  is  thereby  produced  (here,  no  significant  improvement  In  the  location 
of  the  center  of  mass  is  to  be  expected). 

SAGA  also  differs  from  the  other  programs  in  that  more  comprehensive  error 
models  are  employed  for  optical  tracking  systems.  In  addition  SAGA  is  expressly  designed 
to  accommodate  observations  made  by  Geoceivers.  Ranging  error  models  that  have  so  far 
been  incorporated  into  GDAP  and  NAP  are  not  sufficiently  general  to  accommodate 
Geoceiver  observations  (should  the  need  arise,  however,  they  could  readily  be  extended 
to  do  so).  In  the  next  two  sections  we  shall  go  into  the  detailed  development  of  the  optical 
and  electronic  error  models  employed  in  SAGA. 

From  the  foregoing  review  it  can  be  appreciated  that  SAGA  provides  a  powerful 
tool  for  satellite  geodesy.  What  is  not  yet  apparent  is  the  fact  that  in  spite  of  its  flexibility 
SAGA  has  been  designed  to  be  easy  to  set  up  and  use.  This  is  accomplished  in  part  by 
building  into  the  program  selectable  sets  of  standard  options  sufficiently  broad  to  cover  most 
routine  situations  likely  to  be  encountered  in  practice.  Special  situations  can  be  accommodated 
when  required,  but  at  the  expense  of  a  more  extensive  set  up  of  control  parameters. 
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2-°  OPTICAL  OBSERVATIONAL  EQUATIONS 

2. 1  Projective  Equations 

If,  as  in  Figure  1,  X,Y,Z  denote  the  space  coordinates  of  a  point  photographed 
by  a  camera  1c  -  J  at  X‘,  Y^Z*,  it  is  well  known  that  the  image -coordinate.  x,y  are 
ideally  given  by  the  projective  equations 

x  =  x  +c  A<X-XC)  +B(Y-YC)  +C  (Z-Zc) 

D(X-XC)+£(Y-YC)  +  F{Z-ZC) 

^  A'(X-XC)  +  B'  (Y-Y6)  -tC(Z-Zc) 

y  =  yp  +  c  * - - 1 - L 

D  (X-Xc)  +  E  (Y-Y6 )  +  F  (Z-Zc ) 

in  which 


x„/y,,c 


ABC 
A1  B'  C 
D  E  F 


coordinates  of  center  of  projection  is  image  space  (note:  the 
z  axis  of  image  space  is  directed  along  the  camera  axis;  thus 
c  corresponds  to  focal  length). 

orientation  matrix  =  matrix  of  direction  cosines  of  x,y,z  of 
image  space  relative  to  X,  Y,Z  axes  of  object  space 
(specifically,  A,B,C  are  direction  cosines  of  x  relative  to 
X,Y,Z;  A'B'C1  of  y  relative  to  X,Y,Z  and  D,E,F  of  z 
relative  to  X,Y,Z). 


The  orientation  matrix  can  be  expressed  uniquely  in  terms  of  three  angles.  If  the  x,y,z 
axes  of  image  space  and  the  X,  Y,Z  axes  of  object  space  are  related  by  the  angles 
a,  a x  indicated  in  Figure  2,  the  orientation  matrix  can  be  shown  to  assume  the  fonn: 

sin  a  coax  -  cos  a  ainw  tin*  ooawain<c~ 
•in  a  sJ.ax-  cot  a  ainw  cock  ooaweotx 
SOI  a  cot  w  |tg  v 


(2) 


'a  b  c' 

r 

A'  8*  C* 

■ 

D  E  P 

•  m 

oota  •in*  -  tin  a  tin  w  ooax* 
•in  a  aoa<» 


By  dividing  the  numerator  and  denominator  of  the  ratio  on  the  right  hand  of  (1)  by: 
(3)  R  «  [  (X-Xc  Y>  +  (Y-YC)3  +  (Z-Z=  )a 

one  can  express  the  projective  equations  in  terms  of  direction  cosines: 
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(4) 


AX  +  Bp  +  Cv 
x  =  x,  +  c  DX  +  Ep  +717 

A1  X+B'p+C'y 

y  =  yn  +  c  - — ■ 

p  DX  +  Ep+Fv 

wherein 

(5)  X  =  (X-Xc)/R,  /i»(Y-Ye)A  t/=  (Z-Zc)/R. 

Equations  (1)  and  (4)  can  be  put  into  the  alternate  form: 

X-X°  _  A(x-xp)  +  A ‘(y-y,)  +  Dc 
Z»ZC  C(x-xp)  +  C'(y-y„)  +  Fc 

Y-Y*  B(x-xp)+B'(y-yp)+Ec 
Z-Zc  "  C(x-xp)  +C  (y-yp)  +  Fc  ' 

The  direction  cosines  X,p,y  can  be  expressed  also  as: 

X  =  sin  a*  cos  co* 

(7)  p  =  cos  a*  cos  to* 
y  =  sin  to* 

in  which  a*,  to*  are  measured  in  the  same  sense  as  the  a,  to  that  define  the  direction 
of  the  camera  axis  (Fig.  2)- 

Once  the  projective  parameters  a,  t o,x,  xp,yp,c  (and  possibly  others,  such  as 
coefficients  of  distortion)  have  boen  determined  from  a  plate  reduction  based  on  measured 
plate  coordinates  of  selected  stars,  the  plate  coordinates  of  satellite  images  can  be 
employed  in  equations  (6)  and  (7)  to  establish  their  directions  a*,  co*.  If  X,  Y,  Z  are 
suitably  defined,  a*,  co*  will  be  equivalent  to  Greenwich  hour  angle  and  declination, 
which,  in  turn,  can  be  converted  into  right  ascension  and  declination  if  the  time  of 
the  observation  is  known. 

It  has  turned  out  that  optical  observations  published  by  the  various  data 
gathering  organizations  consist  of  the  derived  quantities  right  ascension  end  declination 
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(accompanied  by  time)  instead  of  the  original  observations,  namely,  the  measured 
plate  coordinates.  Thus  the  uncritical  user  is  likely  to  regard  right  ascension  and 
declination  as  directly  observed  quantities  rather  than  as  derived  quantities.  While 
errors  in  plate  coordinates  'remain  uncorrelated  for  all  regions  of  the  celestial  sphere, 
errors  in  right  ascension  and  declination  become  highly  correlated  in  polar  regions. 

Since  most  organizations  do  not  accompany  their  published  right  ascensions  and 
declinations  with  covariance  matrices,  such  correlation  is  not  generally  taken 
properly  into  account.  Consider,  for  example,  the  extreme  case  of  a  plate  centered 
at  the  pole.  If  plate  measuring  accuracies  are  equal  in  x  and  y  (cr„  =  <Jf  -  a),  it  can 
be  shown  that  the  standard  deviations  of  the  derived  right  ascension  and  declination 
(a,  5)  are  given  by: 

Oa  -  (a  tan  6)/c 

(8) 

05  =  (a  sin  a  6)/e 

and  the  correlation  between  a,  6  Is  given  by: 

(9)  -  2  sin  aeos  a. 

Thus  correlations  between  a,  6  can  range  between  -1  to  +1  for  points  on  the  same  plate. 

To  those  versed  in  analytical  photogrammetry,  there  is  good  reason  to  prefer 
plate  coordinates  over  derived  angles.  Aside  from  the  matter  of  correlation,  the 
projective  equations  (1)  relating  x,y  and  X,Y,  Z  are  actually  simpler  than  the  relations 
between  a,  6  ana'  X,Y,  Z.  However,  the  overriding  reason  for  preferring  the  projective 
equations  has  to  do  with  error  modeling.  Systematic  errors  in  optical  directions  are  in 
large  part  attributable  t  rors  in  the  projective  parameters  produced  by  the  plate 
re.  :tion.  Especially  significant,  in  many  instances,  is  the  angular  instability  of 
the  camera  throughout  the  data  gathering  period.  As  will  shortly  be  demonstrated,  a 
physically  meaningful  optical  error  model  can  be  expressed  in  an  especially  compact 
form  when  the  projective  equations  provide  the  observational  equations  for  the  reduction 
In  particular,  we  shall  show  that  four  error  coefficients  can  account  for  a  total  of  eight 
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distinct  sources  of  systematic  error. 

8ecause  plate  coordinates  and  associated  projective  parameters  are  not  generally 
available,  one  may  resort  to  what  may  be  termed  the  'dummy  camera  method*  to 
reconstruct  from  the  given  angles  sets  of  plate  coordinates  that  are  approximately 
equivalent  to  those  actually  measured.  The  dummy  camera  method  Involves  the 
following  steps: 

1)  a  focal  length  c  is  adopted  that  approximates  the  focal  length  of  the 
camera  actually  used; 

2)  a  central  ray  is  selected  to  define  the  direction  of  the  camera  axis 
(oi,  u>)  in  equation  (2); 

3)  with  the  swing  angle  x  provisionally  set  equal  to  zero,  the  orientation 
matrix  of  the  dummy  camera  is  evaluated  from  (2); 

4)  the  given  angles  for  each  ray  are  converted  into  direction  cosines 
X,^,u  by  means  of  such  relations  as  (7)  (typically  a*  corresponds  to 
Greenwich  hour  angle  computed  from  sidereal  time  and  right  ascension  and 
in*  corresponds  to  declination); 

5)  with  Xp  and  yp  set  to  zero,  with  c  equal  to  the  value  adopted  in  step  (1) 
and  with  the  orientation  matrix  computed  in  step  (3)r  the  direction  cosines 

X,  H,V  are  substituted  into  eqs.  (4)  to  generate  equivalent  plate  coordinates 
x,y. 

The  dummy  plate  coordinates  thus  generated  together  with  the  adopted  projective 
parameters  of  the  dummy  camera  provide  artifical  observations  having  errors  equivalent, 
for  all  practical  purposes,  to  the  errors  in  the  original  observations.  For  reasons  shortly 
to  be  made  clear,  one  extra  step  in  fie  dummy  camera  projection  is  desirable.  This  is 
to  redefine  the  swing  angle  x  (provisionally  assigned  a  value  of  zero  In  step, (3))  so  that 
the  x  axis  coincides  approximately  with  the  trace  of  the  satellite,  in  this  regard  we 
would  note  that  when  the  original  observations  are  uncorrelated  and  are  of  the  same 
accuracy  (ax  =  -  a),  so  also  are  transformed  values  x',y'  defined  by: 

x*  =  x  cos  x  -  y  sin  x 
y'  =  y  sin  x  +  x  cos  x  . 
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This  means  that  the  x  axis  of  the  dummy  camera  can  be  arbitrarily  directed  without 
significantly  altering  the  error  structure  of  the  dummy  observations*  Accordingly, 
directing  the  x  axis  of  the  dummy  camera  along  the  traco  of  the  satellite  is  altogether 
acceptable,  even  though  this  may  not  necessarily  correspond  to  the  direction  of  the 
original  x  axis. 


2.2  Optical  Error  Model 

The  x,  y  coordinates  reconstructed  by  the  dummy  camera  method  may  be 
represented  as: 


(10) 
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x  =  x  +  c  - 

p  q 

n 

y  -  yp  +  c  - 


in  which 


m 

A  B  C“ 

"XT 

01) 

n 

_q_ 

s 

A1  B'  C 

D  E  F 

V 

The  systematic  errors  in  x  and  y  attributable  to  systematic  errors  in  projective  parameters 
may  be  represented  as: 


6x  -  6x_  +  —  6c+  —  6m  «*  c  6q 

1  q  q  q3 

(12) 

6y  =  +  ^  5n-c^  6q. 


The  errors  6m,  5n,  6q  arise  from  errors  in  the  orientation  matrix.  Let  Aa,Au>,  Ax  denote 
three  infinitesimal  rotations  that  serve  to  correct  the  orientation  matrix.  Then  if  m^n^q' 


denote  the  correct  values  of  m,n, 

q,  one  can  write: 
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Here,  if  will  be  noted  that  the  matrix  In  the  quantities  Aa,  6  1 0,  6k  qualifies  as  an 
orrhogonal  matrix  if  terms  of  second  order  are  neglected  when  the  matrix  is  multiplied  by 
Its  transpose .  Therefore,  it  is  a  rotation  matrix.  However,  6a,  Aw,  Ax  do  not  consist 
oc  .irect  additive  corrections  to  the  a,  40,  X  implicit  in  the  original  orientation  matrix. 
Physically,  6a  and  Aco  are  components  of  rotation  of  the  camera  axis  in  the  xz  and  yz. 
planes  of  Image  space,  and  Ax  is  a  component  of  rotation  about  the  camera  axis.  The 
errors  in  m,n,q  attributable  to  errors  in  the  orientation  matrix  are  given  by: 
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Before  substituting  these  results  into  (12),  we  shall  express  Aa,  Aw,  Ax  as: 

6a  -  6a  +  T  6a 

(15)  Aw  =  6w  +  r  6w  , 

A  x  =  6x  +  r  6x 

in  which  r  denotes  the  time  of  the  observation  relative  to  the  time  of  the  central  ray 
defining  the  direction  of  the  axis  of  the  dummy  camera.  By  virtue  of  (15),  we  adopt 
the  assumption  that  the  orientation  of  the  camera  is  not  necessarily  strictly  stationary 
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but  may  be  changing  infinitesimally  with  time.  If  we  now  substitute  (15)  into  (14) 
and  then  substitute  this  result  into  (12)  we  shall  arrive  at  the  expressions: 

6x  =  6xp  +  “  5c  -  c  5a  +  6co  +  y  6x 

-c  ^l+4j-jr6a  +  ^ r6w  +  yr  6x 

(16) 

6y  ■  6yp  +  ^  5c  +  ^  6a  +  c  ^1  +  6to  -  x  6x 

+  S^rfia  +  c  ^1  +  t5o >  -  xr  6x  • 

In  the  reduction  leading  to  this  result  we  employed  the  relations  x  =x(m/q),  y  =c(n/q) 
which  follow  from  the  consideration  that  xp  =yp  =0  in  dummy  camera  projection. 

As  it  stands,  the  error  model  (16)  involves  a  total  of  nine  parameters.  However, 
the  number  can  be  reduced  to  a  total  of  four  essential  parameters  by  certain  considerations. 
First  we  note  that  for  cameras  of  long  focal  length  such  as  the  MOTS  and  PC-1000,  x  and 
y  are  less  than  one  tenth  as  great  as  c.  Accordingly,  terms  xa/c3  and  y^/c3  may  be  set 
equal  to  zero  without  significant  effect.  We  now  recall  the  fact  that  the  swing  angle  x 
is  chosen  in  the  adopted  method  of  dummy  camera  projection  so  that  the  x  axis  coincides 
very  nearly  with  the  photographic  trace  of  the  satellite  (which  in  turn  typically  departs 
from  linearity  by  only  a  few  hundred  microns  at  most).  Thus  for  all  points  on  the  trace 
» 0.  Moreover,  the  x  coordinates  of  points  along  the  trace  can  be  represented 
approximately  by  the  relation  x=xr  where  x  denotes  the  mean  rate  of  change  of  x  over 
the  plate.  If  in  line  with  these  cons ide rati ons  we  make  the  following  set  of  sub- 
stitutions  into  (16): 

(17)  xs/c3  =  yVc3  =0,  y  =  0,  xexr 

/ 

we  shall  obtain  the  result: 
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(18a)  fix  =  fix-cfia  +  —  6c  -  cr  6a 
'  c 

(18b)  fiy  =  6y?  -  cfico  -  xr  6  k  -  ct  6  w  -  xtn  fix. 


We  see  in  (18a)  that  the  coefficients  of  6xj,and  fia  are  constant  multiples  of  each  other; 

the  same  is  true  of  the  coefficients  6c  and  6a ,  This  means  that  6a  alone  is  sufficient  to 

account  for  the  combined  effects  of  infinitesimal  changes  in  a  and  xp .  Likewise,  6c 

alone  is  sufficient  to  account  for  the  combined  effects  of  an  infinitesimal  change  in  scale 

(or  focal  length)  and  an  infinitesimal  rate  of  change  in  the  6a  component  of  rotation. 

Similarly  in  (18b)  we  find  that  6yt  and  ficuare  perfectly  coupled,  as  also  are  fix  and 

6 cib •  Thus,  the  rotations  fiui  and  fixserve  also  to  account  for  6yp  and  6w  respectively. 

Further  simplification  can  be  achieved  f  jm  consideration  of  the  fact  that  with  cameras 

/ 

having  a  focal  length  that  is  many  times  larger  than  the  plate  format,  the  term  in  6xis 
likuiy  to  be  relatively  insignificant  in  comparison  with  the  terms  in  6a  and  6cu.  For 
cameras  such  as  MOTS  and  PC-1000  the  coefficients  cr  of  6a  and  6w  in  (18a)  and 
(1 8b)  are  about  ten  times  larger  than  the  maximum  value  of  the  coefficient  xtn  =  xr  of 
fix.  This  means  that  fix  must  be  about  ten  times  greater  than  6a  and  6co  in  order  to 
induce  a  comparable  error.  In  a  study  of  camera  stability  reported  by  Brown  (1969)  the 
maximum  values  of  6a,  6 ii,  and  fix  for  a  PC-1000  were  found  to  be  about  O'.'OI/sec  for 
6a  and  6toand  about  QVQ2/sec  for  6x.  Although  fix  did  become  about  twice  as 
great  as  6a,  fiaj  its  net  effect  was  only  one  fifth  as  great  inasmuch  as  (max  x|  w  0.1  c. 
In  view  of  such  consider. rions,  we  regard  carrying  5x  in  the  error  model  to  be  generally 
of  dubious  value  and  accordingly  have  dropped  It  in  further  treatment  of  the  model . 


By  virtue  of  the  findings  of  the  previous  paragraph,  we  may  drop  from  the 
general  error  model  (16)  the  terms  in  6xp  6yp,  6a,  ficu  and  6x.  This  leaves  a  four 
parameter  model  of  the  form: 


_ 


This  compact  model  is  sufficient  to  account  not  only  for  biases  in  six  projective 
parameters  but  also  for  my  uniform  drift  of  the  camera  axis  throughout  the  exposure. 

One  must  not  lose  sight  of  the  fact  that  this  result  does  depend  in  part  on  dummy  camera 
projection  that  places  the  trace  of  the  satellite  through  the  plate  center  and  approximately 
along  the  x  axis  of  the  plate. 

When  optical  systems  are  employed  to  record  a  flashing  light,  synchronization 
of  all  observations  is  automatic  and  the  problem  of  interstation  timing  bias  does  not  arise. 
However,  when  shutters  are  employed  to  chop  the  traces  of  sun  illuminated  passive 
satellites,  the  possibility  does  arise  that  local  clocks  may  be  inadequately  synchronized. 

In  this  case  the  error  model  (19)  must  be  augmented  by  terms  of  the  form: 


6xt  =  x  6t 

(20) 

6yt  =  y  6t, 

where  fit  represents  the  interstation  timing  bias.  In  c  short  arc  tracking  network  fit  can 
and  should  be  forced  equal  to  zero  for  one  arbitrarily  selected  station  in  the  network. 

The  biases  fit  for  the  remaining  stations  are  subject  to  a  priori  constraints  appropriate  to 
the  timing  system  employed. 

In  cases  where  optical  and  electronic  systems  both  track  a  satellite  carrying  a 
flashing  light,interstation  timing  bias  is  accommodated  in  SAGA  by  treating  the  timing  of 
the  optical  system  es  unbiased  and  the  timing  of  the  electronic  systems  as  biased  relative 
to  the  optical  system. 

A  common  and  desirable  practice  in  optical  tracking  is  to  reorient  each  camera  , 
one  or  more  times  during  the  course  of  a  pass  in  order  to  obtain  extended  coverage  from 
a  given  station,  MOTS  cameras  occasionally  obtain  as  many  as  four  piates  on  a  given 
pass  and  two  or  three  plates  are  common.  Under  such  circumstances  it  becomes 
necessary  to  reinitialize  the  error  model  for  each  plate  (except  for  intersfation  timing 
bias  which  would  be  common  to  all  plates  at  a  given  station  for  a  given  pass).  This 
means  that  if  a  particular  station  were  to  acquire  four  plates  on  a  given  pass,  one 
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would  have  to  determine  an  independent  set  of  coefficients  6  a,  6a>,  6x,  6iC  for  each 
plate  and,  where  applicable,  e  single  interstation  timing  bias  6t  for  all  piates.  As  a 
consequence,  an  optical  station  can  require  the  exercise  of  as  many  as  seventeen  error 
coefficients  for  a  single  pass.  Such  a  capability  is  provided  in  SAGA.  Let  us  consider 
what  this  implies  in  view  of  the  fact  that  SAGA  is  designed  to  accommodate  as  many  as 
fifteen  stations  on  a  given  pass.  The  most  extreme  situation  would  be  one  in  which  all 
fifteen  stations  are  employed  in  a  chopping  mode  and  each  station  successfully  acquires 
four  plates.  The  number  of  error  parameters  to  be  recovered  on  a  single  pass  would  rhen 
amount  to  15x  17 - 1  =  254  (the  timing  bias  at  one  station  is  constrained  to  zero).  Such 
c  reduction  becomes  practical  only  by  virtue  of  the  use  of  second  order  partitioned 
regression  as  is  discussed  in  Section  4. 

2.3  A  Priori  Constraints 

By  virtue  of  the  stellar  control  employed  in  plate  reductions,  systematic  errors 
in  (  "t'ly  determined  directions  are  sharply  bounded.  The  error  budget  for  a  PC-1 000 
reduction  provided  in  Table  2  is  taken  from  Brown,  Bush,  Sibol  (1963).  For  current 
validity  the  budget  need  be  changed  in  only  a  few  respects.  The  use  of  the  SAO  star 
catalog  in  place  of  the  Boss  catalog  would  about  halve  the  contribution  of  item  A3. 
Tangential  or  lens  decentering  distortion  is  now  roulinely  calibrated  and  removed 
according  to  methods  developed  in  Brown  (1964),  (1966).  As  a  result  items  A6  and 
B6  of  the  error  budget  can  be  reduced  to  about  one  third  their  former  values.  The  most 
significant  change  to  the  budget  affects  item  All  which  is  concerned  with  camera 
stability.  The  budget  calls  for  rejection  of  the  plate  if  comparison  between  pre  and 
post  orientations  indicates  the  presence  of  camera  instability  equivalent  to  more  than 
one  third  the  net  rms  error  in  the  p!ate  coordinates.  This  recommendation  has  been 
found  to  be  too  stringent  tc  be  followed  in  general  practice.  Instead,  Instability  is 
tolerated  to  the  point  where  its  effects  on  direction  are  comparable  with  those  random 
errors.  In  effect,  this  means  that  some  PC-1000  plates  are  accepted  even  though  a  change 
In  orientation  of  as  much  as  two  seconds  of  arc  exists  between  pre  and  post  calibrations. 
When  such  a  change  is  continuous  (as  opposed  to  a  suddan  disturbance),  its  effect  on 
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TABLE  2  ERROR  BUDGET  FOR  PC-1000  FOR  POINTS  OUTSIDE  ATMOSPHERE 


Z  Z 

i  2  ^ 

(J  X 
UJ  2  w 

E*  3 

Q  Q  u_ 


RMS  CONTRIBUTION  IN  MICRONS  UNDER. 

(o) 

(b) 

(c) 

ERROR  SOURCE 

Favorable 

Normal 

Unfavorable 

Conditions 

Conditions 

Conditions 

1.  Random  setting  error  (average  of  2  settings). 

1.0 

1.5 

2.0 

2.  Emulsion  instability. 

1.0 

1.5 

2.5 

3.  Low  frequency  atmospheric  shimmer. 

0.5 

1.0 

3.0 

4.  Star  catalog  error  (Boss). 

2.0 

3.0 

4.0 

5.  Residuai  •  adi al  distortion. 

0.5 

1.0 

2.0 

6.  Residual  tangential  distortion. 

1.0 

2.0 

4.0 

7.  Flatness  of  surface  at  emulsion. 

0.0 

0.0 

0.5 

8,  Residual  differential  refraction 

0.0 

0,5 

1.0 

9.  Residual  comparator  errors. 

0.5 

1.0 

1.5 

10.  Timing  errors  (WWV). 

0.5 

1.0 

2.0 

11.  Camera  instability  (below  threshold  of 

routine  detectability). 

1.0 

1.5 

2.0 

POOLED  RSS  TOTALS:  . 

3. Op 

4.9p 

8.2  p 

1 .  Random  setting  error  (average  of  2  settings). 

1.0 

1.5 

2.0 

2.  Emulsion  instability. 

1.0 

1.5 

2.5 

3.  High  frequency  atmospheric  shimmer. 

1.0 

2.5 

5.0 

4.  Residual  error  In  calibrated  orientation. 

0.5 

1.0 

1.5 

5.  Residual  radial  distortion. 

1.0 

1.5 

2.5 

6.  Residual  tangential  distortion. 

1.0 

2.0 

4.0 

7.  Flatness  of  surface  of  emulsion. 

0.0 

0.0 

0.5 

8.  Residual  parallactic  refraction. 

0.5 

1.0 

1.5 

9.  Residual  comparator  errors. 

0.5 

1.0 

1.5 

POOLED  RSS  TOTALS: 

2.4  p 

4.5  p 

8.0  p 

GENERAL  QUALIFICATIONS: 

o.  Calibration  is  assumed  to  involve  at  least  40  stellar  images  compoctly  distributed  about  flashing  light  trace 
and  divided  approximately  equally  between  pre-  and  postcalibrations. 

•  b.  Elevation  angle  of  camera  is  taken  as  30°  and  altitude  of  flashes  as  400nm. 

c.  Photopiocessing  procedure  recommended  by  Gallnow  and  Hageman  (Astronomical  Journal,  pp.  399-404,  Vol« 
61,  Nov.  1956)  is  assumed  employed  in  order  to  minimize  emulsion  instability;  for  same  reason,  points 
within  one  centimeter  of  edge  of  plate  are  assumed  npt  to  be  measured. 

d.  Atmospheric  shimmer  is  taken  to  be  that  characteristic  of  maritime  subtropical  atmosphere  at  30*  elevation 
angle  with  PC-1000  employed  at  full  (200mm)  aperture. 


e.  Timing  errors  (WWV)  are  taken  as  5,  10,  and  20  millisecondsfor  coses  (o),  (b),  (c),  respectively. 

f.  Comparator  is  assumed  to  be  calibrated  and  properly  operated. 

g.  Plates  ore  assumed  rejected  if  comparison  between  individual  pre  ond  post  orientations  indicate 

presence  of  cam'  jility  equivalent  to  more  than  one  third  the  net  rms  error  in  the  plate  coordinates. 


turn- 
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satellite  directions  is  likely  to  be  less  than  one  second  of  arc  because  of  the  relatively 
short  time  interval  spanned  by  the  satellite  observations. 

In  view  of  current  practice,  we  would  suggest  that  a  prior!  constraints  for 
PC-1000  and  MOTS  optical  error  coefficients  be  selected  from  one  of  the  following 
three  schedules: 


Schedule 

Criterion 

°ol  aa>  <*x  °c 

1  (favorable) 

a,  <3|i 

0.'5  C'.'5  r.O  10m 

2  (normal) 

3<<7t<6fi 

U'O  T'.O  2'.'0  20m 

3  (unfavorable) 

o,>6n 

1  "5  1*15  31*0  30m 

The  quantity  a,  refers  to  the  ms  error  achieved  in  the  plate  reduction. 

A  primary  advantage  of  the  short  arc  approach  to  satellite  geodesy  over  the 

geometric  approach  is  the  practibility  of  accounting  for  systematic  errors  in  extensive 

networks  through  error  modeling.  On  strongly  observed  arcs,  adjusted  values  of  many 

of  the  error  coefficients  can  constitute  substantial  improvements  over  a  priori  values. 

On  the  other  hand,  some  error  coefficients  may  prove  to  be  intractible,  their 

accuracies  after  adjustment  being  no  better  than  before  adjustment.  This  has  been 

used  as  an  argument  against  the  exercise  of  error  models  in  the  adjustment.  Such  an 

argument  Is  unsound  for  it  is  clearly  important  that  the  effects  of  statistically  bounded 

systematic  errors  be  rigorously  taken  into  account  even  when  such  errors  are  not 

amenable  to  worthwhile  reduction.  This  is  especially  so  inthe  case  of  plates 

containing  a  large  number  of  satellite  Images,  as  when  passive  satellites  are  recorded. 

Here,  one  wou'd  obtain  unduly  optimistic  estimates  of  accuracy  from  error  propagation 

« 

if  one  were  to  ignore  the  possible  existence  of  systematic  error. 


2.4  Special  Corrections 

SAGA  is  designed  to  accept  any  optical  observations  that  are  provided  in  the 
GEOS  format  of  the  NASA  data  bank.  Unfortunately,  toere  is  no  uniform  standard 


with  regard  to  certain  corrections  that  is  followed  by  oil  agencies  producing  optical 
data.  In  particular,  the  following  corrections  may  or  may  not  have  been  applied  by  a 
given  agency: 

1)  polar  motion; 

2)  conversion  of  times  to  UT1 ; 

3)  parallactic  refraction; 

4)  phase  angle  correction  (chopping  of  passive  satellites); 

5)  propagation  delay  (chopping  of  passive  satellites). 

Because  of  the  lack  of  homogeneity  in  the  application  of  such  corrections,  SAGA  has  been 
provided  with  a  special  preprocessor  which  serves  to  apply  these  corrections  as  needed. 
Characteristics  of  the  preprocessor  are  discussed  in  Part  II  of  this  report. 


2.5  Orbital  Constraints 


The  orbital  integrator  employed  by  SAGA  is  that  developed  by  Hartwell  (1968), 
(1969).  It  employs  a  power  series  solution  to  the  equations  of  motion  wherein  each 
coefficient  is  formed  in  terms  of  its  predecessors  by  means  of  recursive  algorithms.  The 
version  of  the  integrator  employed  in  SAGA  truncates  the  gravitational  potential  at 
(n,m)  =  (4,4)  inasmuch  as  this  has  proven  to  be  entirely  adequate  in  short  arc  applications 
(Brown  1967).  If  x,y,z  denote  the  geocentric  inertial  coordinates  at  an  arbitrary  time  r 
relative  to  an  adopted  epoch  r=0,  the  power  series  solution  of  the  equations  of  motion 
can  be  represented  as: 
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in  which  all  of  the  coefficients  are  functions  of  the  six  initial  conditions  at  r  -  0 
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(namely:  x0/y0,zo  anc*  gravitational  coefficients.  The  series  is  truncated 

automatically  when  a  prespecified  tolerance  (presently  taken  as  0,001  m)  is  satisfied 
for  the  maximum  value  of  r  to  be  exercised.  If  the  epoch  is  taken  near  midarc,  the 
radius  of  convergence  of  each  expansion  is  suf'iciently  great  to  accommodate  arcs  as 
long  as  one  third  of  a  revolution  for  nearly  circular  orbits. 

The  version  of  the  integrator  employed  in  SAGA  also  geneiates  power  series 
solutions  to  the  variational  equations  relating  errors  in  x,y,z  at  time  r  to  errors  in  the 
adopted  location  of  the  center  of  mass  and  errors  in  the  six  initial  conditions.  If  we 
let  Xqq ,  Ym  ,  Zm  denote  the  earth-fixed  coordinates  of  the  center  of  mast,  the  matrix 
of  partial  derivatives  generated  by  the  integrator  (the  matrizant)  can  be  expressed  as 


(22) 
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in  which  each  is,  in  turn,  a  polynomial: 
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Inertial  coordinates  generated  by  the  orbital  integrator  can  be  referred  to  an 
earth-fixed  framework  by  the  application  of  the  transformation: 
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in  which 
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cos  l (>T  sin  0 

-sin  ipT  cos^t  0 

(25)  R  = 
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0  0  1_ 

°  0  0 

in  a  similar  manner  the  matrisant  can  be  transformed  to  earth  fixed  coordinates  by  the 
operation: 


(26)  $  *  =  R  Ci 

(3,0)  (3,  3)  (3,  a) 


6(X,Y,Z;t) 

_  '  e 

^(^00  >  Yoo  >  ^00  /^o  »yo/  ^o'^yo'^o) 


With  these  results  we  are  now  in  a  position  to  develop  the  optical  observational  equations 
for  the  short  arc  reduction. 


2.6  Optical  Observational  Equations 

If  the  X,  Y,Z  in  the  projective  equations  (1)  are  replaced  by  the  values  computed 
from  (24),  the  equations  become  functionally  of  the  form: 

x  =  f1(X^Y«,ZS)^Y<#,Z*#i^y0,^y<)#*0f-t) 

(27) 

y  —  fa  (XC,YC,  Zc ,  Xqq  ,  Ygg ,  Zqq  ,  x0,  yo'^o'Xo  •  Yo>  »t)« 

If  x°,y°  denote  the  observed  values  of  x,y,  the  edjusted  values  corrected  for  systematic 
error  can  be  expressed  as: 

x  =  x°  +  v,  +  x  vt  +  x  8t  +  fix 

(28) 

y  =  y°  +  vy  +  y  vt  +  y  6t  +  6y 

where  6x,  6y  are  given  by  (19),  In  (28)  v„,vy  denote  residuals  reflecting  random  error 


^  SiHMWWs 
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in  x  and  y,  and  vt  denotes  the  residual  reflecting  random  error  timing.  The  terms 
in  5t  account  for  the  bias  in  timing  (eq.  (20)).  We  now  set  up  the  relations: 

Xc  =  (Xc)°°  +  6XC  Xgo  =  X£5  +  6X00  x0=^°+6x0  i0=x«+6i0 

(29)  YC=(YC)'>=>  +  6YC  Y00=YS+6Y00  y0  =  yS°  +  6y0  yo'tf'+Syo 

Zc  =(Zc)oo+ 6Zc  2<0=2S+6Z0j  ib=3°+fi*o 

in  which  the  superscripts  (oo) denote  approximations  and  the  6's  ore  corresponding 
corrections.  Substituting  these  into  the  right  hand  side  of  (27)  and  linearizing  the 
resulting  expressions  by  Taylors  series,  we  obtain: 


x  =  x00  + 


3x 


S(XC ,  Yc  ,ZC/XCD ,  Yqq  ,  Z ,Xq  ,  yo/Ze/Xo/yafZo) 


(30) 


y  =  y00 
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d(Xc,Yc  ,Zc,X0o,  Y00/Z00,x0,y0/z0,x0,y0,z0) 


(6XC ,  AY®,  GZC , . . fiy0 


(8XC ,  6 Y6, 5Ze ,  •  •  •  #  8y0 


in  which 


(31) 


x«  *  fxWX6)00,^)00,...,*?1) 
y00  =  fa((Xc)00,  (Yc)°° . z°°). 


To  arrive  at  explicit  expressions  for  the  elements  of  the  linearized  observational  .equations, 
we  let  X00,  Y10 ,ZM ,X°° ,Y30,Z00  denote  the  components  of  position  and  velocity  for  the 
time  T  of  the  observation  as  computed  from  equations  (21)  and  (24)  in  which  the  given 
approximations  in  (29)  are  employed  in  the  integration.  Then  we  define  the  auxiliary: 
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where  the  orientation  matrix  has  been  developed  by  dummy  camera  projection.  The 
values  of  the  plate  coordinates  computed  from  (31)  then  become: 
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where  c  denotes  the  focal  length  adopted  in  dummy  camera  projection.  The  partial 
derivatives  of  the  plate  coordinates  with  respect  to  station  coordinates  is  given  by: 


(34) 

(a, a) 


8fc,y)  _  _-c  1  '°  **°°/c 

d(Xc,Y«,Zc)  <f°  0  1  -y°°/c 


In  terms  of  this  the  partial  derivatives  of  the  plate  coordinates  with  respect  to  center  of 
mass  and  orbital  state  vector  are  given  by: 

(35)  B<a>  =  $ 

(3/8 )  ^(^oo#y»r  (a, a)  (3  ,9) 


in  which  $  is  given  by  (26).  The  time  derivative:'  of  the  plate  coordinates  required  in 
(28)  can  be  computed  from: 


If  we  partition  as: 
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the  linearized  observation  equations  can  be  put  into  the  form: 

(38)  A  v  -  BW  fid)  -  Bd‘)  fid*)  -  bM  fiW  -  bW  5(3)  .  B«)  6(*)  «  ( 

(a,3)(3/i)  (a,s)(s,i)  (3/3)(3/i)  (a,e)(s,  i)  (a,i)(x,i)  (a,4)(*,i)  (a,i) 
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in  which 
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At  this  point  we  shall  recognize  that  as  many  as  four  plates  may  be  acquired  at  a  given 
station  for  a  given  pass.  Accordingly,  there  may  be  as  many  as  four  sets  of  error 
coefficients.  Letting  p  denote  the  pth  plate  (max  p=4)  recorded  and  introducing  the 
subscript  j  to  denote  the  j  th  point  observed  by  the  station,  we  may  express  the  pair  of 
linearized  observational  equations  generated  j  th  point.  If  observed  on  plate  p,  as; 

(40)  AjVj  +  B}  6  ® 
in  which 
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covariance  matrix  of  the  random  errors  in  plate  coordinates  and  timing  for  the  j  th 
point  from  the  given  station  is  denoted  by: 


(44)  A, 

(3/3) 


the  system  of  normal  equations  generated  by  the  point  can  be  expressed  as: 


(45)  Nj  6  =  Cj 
in  which 

N  j  «  B;<AJA1AJTflBj 

<46) 

c)  ~  Cj . 

It  is  to  be  noted  that  since: 

„3  ,  *3  3 

a>J  +  “l<rT, 

°  V^J 

the  results  of  the  multiplication  by  (A {  A}  Aj)  can  also  be  effected  by  treating  this 
matrix  as  a  unit  matrix  in  (44)  after  modifying  B }  and  e }  by  dividing  their  first  and 
second  rows,  respectively,  by  (cr^  +x*  )^  and  (cr^  +  y °  cr^)^  . 

The  system  normal  equations  generated  by  all  points  from  the  given  station  and 
pass  is  simply: 

(48)  N  6  =  c 

where 

'N  =  LN. 

(48a) 

c  =  2  c,  . 

■  t 
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We  now  introduce  subscripts  i  ond  k  to  denote  the  i  th  station  and  kth  pass  respectively. 
The  normal  equations  (48)  may  then  be  written  in  more  detailed  partitioned  form  as: 


We  shall  find  it  convenient  to  proceed  formally  as  if  all  m  stations  in  the  tracking  network 
were  to  observe  all  passes  (presently,  this  assumption  will  be  dropped)  If  we  then  assume 
that  (49)  has  been  evaluated  for  all  m  stations  (i  =  1,2, . .  .,m),  and  merge  the  resulting 
individual  sets  of  normal  equations  into  a  common  system  by  the  process  of  zero  augmentation 
developed  in  Brown,  Trotter  (1967),  we  shall  obtain  the  system  (50)  indicated  on  the  next 
page.  If  a  particular  station  does  not  participate  in  the  tracking  of  the  kth  pass,  equation 
(50)  should  be  modified  by  (a)  replacing  all  elements  corresponding  to  the  station  in  the 

A 

U  and  c  portions  of  the  normal  equations  by  zeroes,  and  (b)  deleting  from  the  remainder 
of  the  normal  equations  the  elements  corresponding  to  the  station. 

By  adopting  tho  partitioning  indicated  by  the  broken  lines  in  (50),  we  can 
represent  the  system  of  normal  equations  for  the  kth  pass  in  the  more  compact  form: 


We  next  assume  that  (51)  has  been  evaluated  for  all  n  observed  passes  (k  =1,2, . .  .,n) 
and  again  resort  to  the  process  of  zero  augmentation  to  merge  all  such  individual  sets  of 
normal  equations  into  a  common  sysiem.  This  leads  to  the  system  indicated  in  equation 

A  •  *• 

(52)  below.  We  shall  assume  that  weight  matrices  W,  Wk ,  W*  reflecting  the  a  priori 
constraints  to  be  exercised  in  the  adjustment  have  been  absorbed  into  the  appropriate 
diagonal  blocks  of  (52).  Then  (52)  represents  the  final  system  of  normal  equations. 
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Having  fonned  the  system  of  normal  equations  for  the  adjustment,  we  must  now 
address  the  problem  of  solving  the  system,  expecially  in  view  of  the  consideration  that 
it  can  grow  to  huge  dimensions.  We  do  this  in  Section  4  which  is  devoted  to  the 
theory  of  partitioned  regression.  We  need  only  point  out  here  that  the  system  of  normal 
equations  (52)  is  precisely  of  the  same  form  as  the  second  order  partitioned  system 
indicated  in  equation  (66)  of  Section  4.  Accordingly,  by  applying  the  algorithms  of 
second  order  partitioned  regression  as  developed  in  Section  4  to  the  present  problem, 
one  can  execute  a  practical  solution  of  the  normal  equations  no  matter  how  many  arcs 
are  processed  or  how  many  error  parameters  are  introduced  by  each  arc. 


2.8  Analysis  of  Residuals 


After  the  normal  equations  have  been  solved  by  the  algorithm  for  second  order 
partitioned  regression,  the  optical  residuals  can  be  computed  from: 


-  AiA)(Ai\i  A’f1^  -Bj  6). 


These  residuals  are  employed  in  SAGA,  along  with  residuals  corresponding  to  other 
observed  quantities,  to  compute: 

s  =  I  (quadratic  form  of  all  residuals) 
l  (degrees  of  freedom) 


The  observational  equations  are  then  relinearized  about  the  new  approximations  to  the 
parameters  and  the  solution  is  iterated,  leading  to  fresh  residuals  and  a  new  value  for  s* 
This  process  is  repeated  until  successive  values  of  s  differ  by  less  than  a  preset  criterion 
or  the  maximum  allowable  number  of  iterations  have  been  executed.  Computational 
details  are  given  in  Section  6. 

Upon  convergence,  the  final  observational  residuals  from  designated  channels 
can,  on  option,  be  subjected  to  an  autoregressive  analysis.  When  this  option  is 
exercised,  the  entire  solution  is  repeated  with  serial  correlation  being  duly  considered 
in  accordance  with  the  process  of  autoregressive  feedback  developed  in  Section  5. 


2.9  Master  Survey  File 


Before  the  adjustment  of  satellite  observations  can  be  undertaken,  it  is 
necessary  to  set  up  a  Master  Survey  File  for  the  entire  tracking  net  ultimately  to  be 
adjusted.  This  file  contains  the  initial  coordinates  of  all  stations  and  the  covariance 
matrix  of  the  coordinates. 

When  initial  coordinates  are  expressed  as  geographical  coordinates,  a 
preliminary  program  (SET-UP)  transforms  these  coordinates  into  Geocentric  Cartesian 
coordinates  and  computes  the  associated  covariance  matrix.  In  order  to  permit  the 
investigator  considerable  flexibility  in  imposing  interstation  constraints,  SET-UP 
also  permits  the  directional  components  and  the  length  of  the  vector  joining  arbitrary 
pairs  of  stations  to  be  constrained  to  any  desired  degree.  In  addition,  it  permits  the 
introduction  of  any  number  of  linear  constraints  between  arbitrary  pairs  of  stations. 

Thus  if  (X“,  Y^,Z“),  (X^°,  Y^°  ,Z°°)  denote  the  initial  coordinates  of  stations  pand 
q,  the  program  admits  constraints  of  the  form: 

(54)  a1X“  +  aaY«+a.Z“  +  01Xf  +  0a  Y?  +  0a  Z?  =  Up4 

in  which  Up,  is  considered  to  be  an  observation  of  prespecified  variance  CTp4  and  the 
a's  and  /3's  are  prespecified  coefficients.  Such  constraints  serve  a  variety  of  functions, 
including  (a)  holding  selected  stations  fixed  relative  to  a  designated  station  (datum 
constraints),  (b)  defining  directions  of  coordinate  axes  (for  adjustments  limited  to  ranging 
observations),  (c)  imposing  special  relations  between  stations. 

All  specified  interstation  constraints  are  properly  exercised  by  SET-UP  in 

A 

developing  the  a  priori  covariance  matrix  A  of  the  geocentric  coordinates  of  the 
tracking  net.  Details  of  this  process  are  given  in  Part  II.  The  location  of  the  center 
of  mass  with  respect  to  the  adopted  geometric  origin  is  treated  as  if  it  were  the  last 
tracking  station  of  the  net.  By  virtue  of  this  artifice,  the  covariance  matrix  A  serves 
to  accommodate  the  a  priori  covariance  matrix  of  the  coordinates  of  the  center  of  mass. 

The  Master  Survey  File  is  set  up  but  once  for  a  given  tracking  net.  SAGA 
calls  upon  this  file  whenever  it  requires  either  the  a  priori  coordinates  of  a  given  station 

A  — ^ 

or  the  a  priori  weight  matrix  (W  ■  A  )  of  the  entire  vector  of  station  coordinates. 
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3.0  ELECTRONIC  RANGING  OBSERVATIONS 


3.1  The  Geocelver 

In  this  section  we  shall  develop  observational  equations  and  error  models 
appropriate  to  various  electronic  ranging  systems:  lasers,  Secor,  GRARR,  Radars,  and 
Geoceiver.  We  shall  place  particular  emphasis  on  the  Geoceiver  system  both  because 
of  its  potential  future  importance  to  satellite  geodesy  and  because  its  error  model  turns 
out  to  be  sufficiently  broad  to  encompass  those  of  all  other  ranging  systems. 

Characteristics  of  the  Geoceiver  are  discussed  in  Stansell,  et.  al.  (1965). 

Briefly,  the  Geoceiver  is  a  compact,  self-contained,  relatively  inexpensive  (under 
$100K),  manpack,  doppler  tracking  unit  designed  to  track  Transit  satellit  >s  as  well  as 
other  satellites  radiating  on  either  of  the  frequency  pairs:  162/324M  Hz,  150/400MHz. 
Reception  of  dual  frequencies  provides  tl  means  for  correction  of  ionospheric  refraction. 
Like  t'1 .  Tranet  system,  the  Guuceiver  exploits  one  way  doppler.  However,  Tranet 
opiates  on  cycle  counts  in  a  destructive  mode;  that  is,  at  preset  intervals  it  measures 
the  time  required  to  acquire  a  preset  number  of  doppler  cycles,  whereupon,  it  ceases 
counting  until  the  start  of  the  next  interval.  Thus,  continuity  of  cycle  count  is  lost 
and  it  becomes  necessary  to  treat  Tranet  observations  as  being  essentially  a  measure  of 
doppler  frequency  (or,  equivalently,  of  range  rate).  By  contrast,  Geoceiver  operates 
on  cycle  counts  in  a  nondestructive  mode;  that  is,  continuity  of  cycle  count  Is 
preserved  (as  long  as  phase  lock  is  maintained).  It  is  this  fact  that  permits  Geoceiver 
to  be  viewed  as  being  inherently  a  ranging  system,  for  if  the  transmitted  frequency  from 
the  satellite  and  the  reference  frequency  generated  by  Geoceiver  were  perfectly 
matched  and  both  were  perfectly  stable,  the  scaled  cumulative  cycle  count  of  beat 
frequency  would,  except  for  an  unknown  additive  constant  of  integration,  represent  a 
direct  measure  of  slant  range.  In  practice,  this  measure  is  contaminated  by  the  unknown 
offset  between  transmitted  frequency  and  local  reference  frequency  as  well  as  by  any 
drifts  of  the  two  frequencies.  Such  factors  can,  however,  be  taken  into  account  by 
error  modelling. 
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(4) 


X  =  c/f0 
Af  =  fj  -  f  . 
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The  quantity  X  represents  the  wavelength  of  the  traramitted  frequency  and  the  quantity 
represents  the  beat  frequency  generated  by  the  Geoceiver.  If  we  now  assume  that  f0  and 
fa  are  constant  over  the  tracking  interval,  we  can  Integrate  both  sides  of  (3)  between  an 
arbitrary  time  r  =  0  and  a  later  time  Tto  obtain: 

(5)  r-r0  =  XN+X(fJ  -f0)r+Ar. 


in  which 


r  ®  range  at  time  t. 


r0  =  range  at  time  T  =  0, 

N  -  / T  Af  dt, 

0 


Ar. 


JL 

2c 


x 

/  v3  dt  ■  special  relativistic  correction, 
o 


The  integral  defining  N  represents  the  number  of  cycles  of  beat  frequency  accumulated  from 
the  initiation  of  counting  (t*0)  until  time  r . 

In  the  case  of  Geoceiver,  cycle  counts  are  cumulated  over  intervals  of  nominally 
one  minute.  Specifically,  Geoceiver  generates  the  cycle  counts: 


(6)  AN,  =  /Tj  Af  dt 
Tfi 

in  which  t,~  x  and  r,  represent  the  times  of  the  first  positive  cycle  crossings  following 

•  '•  i'.y  ■  -  ’ 

successive  one  minute  marks  T,_x  and  T, ,  as  shown  in  the  accompanying  figure. 
Following  readout,  the  counter  is  reset  to  zero  before  the  next  positive  cycle  crossing 
is  counted  in  the  next  interval  and  continuity  of  cycle  count  is  maintained.  The  cycle 
count  N  appearing  in  (5)  can  be  related  to  the  counts  generated  in  (6)  by  the  relation : 


(7) 


Ts 


N,  =  /  Af  dt  =  f  Af  dt  +  f  Afrit  +  .^+J1 
1  o  o  Ti  Tj.x 


Af  dt. 


=  ANX  +  ANa  +  ...  +  AN, , 
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.  One  might  expect  from  (7)  that  errors  in  successive  N }  are  highly  correlated  by  virtue 
of  being  formed  of  common  ANt .  This  is  not  the  cose,  however,  when  Gooceiver  is 
functioning  properly  (i.e.,  no  cycles  are  dropped  from  or  added  to  each  count)  for 
the  ANj,  being  integers,  may  be  viewed  as  being  free  of  error.  The  quantity  properly 
to  be  regarded  as  subject  to  error  is  T}f  the  time  of  the  end  of  each  counting  interval# 
Errors  in  successive  times  Ti  and  are  not  likely  to  be  strongly  correlated. 

In  practice,  the  frequencies  f0  and  f^are  not  perfectly  known,  nor  are  they 
perfectly  stable.  To  account  for  biases  and  linear  drifts  in  these  frequencies  we  may 


write: 

fo 

=  ^00 

+  fifo 

+  f<»T 

(8> 

<• 

fo 

- 

+  6fo 

+  f«5  r 

in  which 

foo  =  adopted  value  of  transmitted  frequency 

=  adopted  value  of  local  reference  frequency 

6f0  =  bias  in  adopted  value  of  f0  at  t=0 

6fo  =  bias  in  adopted  value  of  f^  at  T- 0 
•  « 

f0,fi  =  drift  rates  of  f0  and  f£,  respectively. 

For  Geoceiver,  the  offset  Af^  =  f®  ”  foo  between  the  adopted  frequencies  can  range 
between  16Kc  and  32  Kc,  depending  on  the  type  of  satellite.  Because  of  this  and  the 
fact  that  doppier  ranges  between  ±10 Kc,  the  beat  frequency  Af  will  never  cross  through 
zero.  This  obviates  the  need  for  distinguishing  between  positive  and  negative  cycler  < 
if  we  substitute  (8)'  into  (3)  and  expand  the  resulting  expression  for  )/f0  into  the  series 


(9) 


& 

^oo 


ia 

*co 


we  shall  obtain  the  result 


(10)  r  »  X0  ^1  -&J  (Af  -  Afco)  +  X0(6f0  -  5f«;>  +  X0(f0  “fi)r  +  \  ~  +  higher  order  term*, 
where 

(11)  X0  =  7-—  =  wavelength  corresponding  to  adopted  frequency  of  transmission. 

.  ‘°° 

Integrating  this  expression/  as  before,  we  get: 

(12)  r  -  r0  ■  X0(l  -T1)  "  Af°°  +  Xo  “  ^o)  T  +  \  xo  tf©~fe>  T“ 

f  Lt%  +  higher  order  terms. 

This  can  be  put  into  the  form: 

(13)  V  '  ■  X0(N  -  Afoo  r)  +  a0  +  alT  +  n3  r>  +  a3r  +  Lrt  +  higher  order  terms 


where: 


a3  “(Gfo/^oo)* 

Substituting: 

(15)  r  =  [(X-Xc)3  +  (Y-Yc)3  +  (Z-Zc)3] 

(16)  r°=  X0(N-Af00  t) 

into  (13)  we  obtain  the  basic  equation 
*  See  equation  (25)  for  explanation. 
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(17)  r°  +a0+alT+aar3  +oar  +  6r.  =  [(X-Xc)3  +  (Y-Yc)a  +  (Z-Zr)3]*. 

The  development  thus  far  has  failed  to  take  into  account  the  refractive  effects 
of  the  ionosphere  and  troposphere.  Nor  has  it  taken  into  account  the  effects  of  (a) 
propagation  delay,  (b)  interstation  timing  bias  and  (c)  general  relativistic  effects. 

Thus  to  render  (17)  more  nearly  correct  we  should  odd  to  r°  the  composite  correction: 

(18)  Ar  =  Arx  +  ArT  +  Arp  +  Arr  +  Ar^ 
where 

Ar£  =  correction  for  ionospheric  refraction 

ArT  =  correction  for  tropospheric  refraction 

Arp  =  correction  for  propagation  delay 

Ar^  =  correction  for  interstation  timing  offset  (or  bias) 

Ara  =  correction  for  effect  of  gravitational  potential  or  frequency  (general 
relativistic  effect). 

The  two  frequency  method  is  used  by  Geoceiver  to  determine  ionospheric 
refraction.  If  AAN,  denotes  the  refraction  cycle  count  (i.e.,  the  cycle  count  of  the 
beat  frequency  between  the  two  received  frequencies),  the  desired  correction  is  given 

by: 

(1?)  Arx  =  -K  X0(AANx  +  AANa  +  ...  +  AAN0) 

where  K  =  1/9  for  the  frequency  pair  162/324  Me  and  K  =  1/9^  for  the  frequency  pair 
150,400Mc.  Alternatively,  the  correction  could  be  effected  by  replacing  each  ANj  J 
in  (7)  by  ANj- K  AANr 

The  correction  for  tropospheric  refraction  can  be  computed  by  the  following 
formula  derived  by  J.  Willmann  (private  communication): 

(20)  ArT  =  -2a  f(E) 

where 
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a  =  (n0  -  1)  H0 

f(E)=  1/  [sin  E  +  (sinaE  +  ^)~] 


in  which 

n0  “  index  of  refraction  at  Nation 
r0  =  radius  of  earth  (meters) 

H0  =  scale  height  of  troposphere  (meters). 


The  value  of  n0  can  be  computed  from  meteorological  measurements  at  the  station  by  means 
of  the  formula- 


where 


P0  =  atmospheric  pressure  (mm/Hg) 
e0  -  water  vapor  pressure  (mm/Hg) 
T0  =  temperature  (deg  Kelvin). 

The  scale  height  H0  can  be  approximated  by: 


(23)  H0  =  29.2  (T0 -30). 

The  correction  for  propagation  delay  serves  to  correct  the  range  of  the  satellite 
to  the  position  occupied  at  time  r  j  as  measured  by  the  satellite  clock  (this  is  equivalent 
to  the  range  at  time  +  —  as  measured  by  the  station  clock): 

(24)  Arp  =  f  r-  . 

v 

When  a  Gecceiver  station  is  considered  as  one  of  several  stations  participating 
in  the  tracking  of  a  pass,  it  becomes  necessary  to  allow  for  the  possibility  of  a  significant 
bias  in  the  correlation  between  the  clocks  at  the  various  stations.  If  6T0  denotes  the 
offset  a*  epoch  of  the  clock  at  a  particular  Geoceiver  station  from  the  master  clock,  the 
offset  at  time  r  is: 
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fif0  ,  6f0 

6r  =  6r0  +- —  r  =  6r0  +  X0 - r 

*00  c 


(for  readout  triggered  by  timing  signals  from  the  satellite  clock),  and  the  correction  to 
be  applied  to  r°  is: 

(25)  ArT  =  r6T  =  rfir0 +X0  6f0  ^t- 

When  the  offset  is  unknown,  6r0  can  be  carried  as  an  unknown  parameter,  subject  to 
appropriate  a  priori  constraints.  The  contribution  of  6f0  can  be  absorbed  into  the  error 
coefficient  ax.  For  readout  triggered  by  the  local  clock,  6f0  in  (25)  should  be  replaced 
by  6f<5 . 

The  final  correction  is  attributable  to  the  effects  on  frequency  of  the  difference 
in  gravitational  potential  between  the  satellite  and  the  receiving  station.  The  correction 

is  given  by: 


(26;  Ar .  -  X 

C  t 


J  (h/ra  )dr 

o  TTTO 


where 


p  =  gravitational  constant*  4x10 14 m3/*®^ 
h  =  altitude  of  satellite 

r0  -  distance  of  station  from  center  of  earth  (same  units  as  h). 

It  is  to  be  noted  that  for  nearly  circular  orbits,  h  is  nearly  constant  over  a  pass.  Consequently, 
for  such  orbits  Ar0  becomes  a  linear  function  of  r,  and  the  term  in  ax  in  (17)  can  absorb  the 
effects  of  the  general  relativistic  correction.  A  similar  remark  applies  to  the  special 
relativistic  correction  Ar,  inasmuch  as  velocity  v  is  nearly  constant  for  nearly  circular  orbits. 

For  greater  flexibility  we  shall  augment  the  Geoceiver  error  mode!  by  terms  to 
account  for  residual  interstation  timing  bias  and  residual  tropospheric  refraction  (we  assume 
that  the  corrections  indicated  above  have  been  applied  but  are  not  wholly  adequate).  Thus 
the  full  error  model  becomes: 

(27)  6r  =  a0  +axr +  aars +o3r +  a4r +aBf(E) 
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where  the  terms  in  a4  and  as  account  for  interstation  timing  bias  and  residual  tropospheric 
refraction,  respectively.  Regarding  residual  tropospheric  refraction  we  would  remark 
that  ray  tracing  exercises  through  actual  atmospheric  profiles  indicate  that  if  the  best 
possible  value  of  ft  were  employed  in  (20),  the  proportional  error  6RT/RT  in  the  tropospheric 
refraction  correction  could  be  held  to  under  one  percent  for  E  >  10°  (or  to  about  0.2  m  at 
E  =  1Q°).  However,  when  ft  is  computed  from  meteorological  measurements  made  at  the 
station,  6RT/RT  is  likely  to  range  from  three  to  five  percent  for  E  >  10°  (and  thus  become 
as  great  as  1 .0m  at  E  =  10°).  Hence,  potential  improvement  by  a  factor  of  from  three  to 
five  is  possible  from  error  modeling. 


i 

f 

ii 
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3 . 3  Linearized  Observational  Equations 

We  have  already  remarked  that  the  quantity  subject  to  error  in  Geoceiver 
observations  consists  actually  of  the  time  associated  with  the  cycle  count  (here,  we 
ignore  momentarily  the  contribut  ion  of  errors  in  corrections  for  ionospheric  refraction). 
Let  us  suppose  that  for  c  given  value  of  r°  in  (17)  the  measured  time  is  r  and  the  correct 
time  is  r  ”  tT  (thus  eT  denotes  the  error  in  r).  Then  we  can  write 

(28)  r°  (t-ct)  =  r°(r)  -  r(r)  +  higner  order  terms. 


Now,  €r  is  attributable  in  part  to  errors  in  phase  measurement  and  in  part  to  errors  in 
quantization.  Thus,  we  may  write; 


(.?9)  €T-4  +  €; 

where 

=  contribution  of  phase  measuring  error 
€"  =  contribution  of  quantization  error. 

If  denotes  the  phase  measuring  error  (expressed  as  a  fraction  of  a  cycle)  giving  rise  to 
*'r,  we  may  write 


(30) 
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tjjTil  ttl  X  I,  ■.»  w*  r-JM  nmtm.  *- 


where  Aiydenotes  the  beat  frequency  at  time  r.  But  from  (10)  Afr*  ~  r  (r),  so  that 

(31)  <■}.=  X0  c^j/i ‘(t). 

Substituting  this  into  (29)  and  substituting  the  resulting  expression  into  (28),  we  get: 

(32)  r°  (r-cT)  «  r°(r)  -  X0  e0  -  r(r)  . 

Letting  vr  and  vT  denote  residuals  corresponding  to  the  errors  X0  and  we  may 
write  the  Geoceiver  observational  equations  in  the  expanded  form: 

(33)  r°  +  vr  +  rvT+  Ar  +a0  +  0^  +  aara  +a3r  +  a4r  +  a6  f(E)  = 

[<X-XC)S  r  (Y  ‘Y0  )3  +  (Z-Zc  )a  . 

If  we  now  replace  X,Y,Z  on  the  rig!  hand  side  of  (33)  with  the  values  computed  from 
the  orbual  integrator  (eqs  (24)  and  (21)  of  Sec.  2),  we  shall  obtain  an  equation  of  the 

functional  form; 


(34)  r  =  f(X=,Y«,Z=,X0),Y„,Z. 
where  r  denotes  the  left  side  of  (33). 

From  this  point  on  we  proceed  precisely  as  we  did  in  Section  3  with  optical 
observations.  We  express  the  unknowns  in  (34)  in  terms  of  the  same  initial  approximations 
and  corrections  as  were  used  in  the  optical  reduction  (eq.  (29),  Sec.  3),  and  thus  expand 
the  resulting  expression  in  Taylor's  series  to  get: 


(35) 


=  r°°  + 


»r 


S(XC ,  Y0,  Zc,  Xjo  ,  Yay  Zqo,  Xo ,  y0 ,  Zy  Xq,  y0 ,  Zo ) 


(fiXS6YS6ZS...,6y0,fii0)T 


in  which 

(36)  r"°  =  f((Xc)00,(Yc)00, .  ..,Zq°, r) . 

If  we  lot  X00,  Y**, Z00,)?00,  Y^Z00  denote  the  earth  fixed  coordinates  of  the  satellite  as 
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determined  from  the  integration  based  on  the  approximate  state  vector,  the  partial 
derivatives  of  r  with  respect  XC,VC,ZC  can  be  expressed  as: 

071  (lf"'5kvT?)=-cx'‘''1 

in  which  < 

X  =  (X00  -  (Xc)°°  )/r°° 

(38)  n  -  (Y°°  -  (Yc  )°°)/r00 
V  =  (Z00  -  (Zc)M)/rM. 

The  partial  derivatives  of  r  with  respect  to  center  of  mass  and  orbital  state  vector  are  then 
given  by: 


(39)  S(3)  =  — 


=  - - r-  -  -  B«  $ 

(1,9)  S(X<d/Yooa  (1, 3) (3,9) 


where  $  is  given  by  (26)  of  Sec,  3.  The  approximate  value  of  r  can  be  computed  from: 


(40)  r00  =  -  B(l)  Y03  . 


If  we  partition  B's  as: 


(4i)  b(s)  =  r~ — — 

(1,9)  .  Yco'Zjo) 


B(*\j  /  Yor  *Ot*0'Yofz<^ 


-1  =  f"  B(3l)  i  B(3b^ 
U  LM)  i(*e)_ 


(he  linearized  observation  equation  can  be  expressed  as: 

(42)  A  v  -  B(lJ  6^  8(s  »)  fi(3  *)-  B<sb)  -  B(3>  6<3>  -  6W  -  < 

(l/3)(3,l)  (l,a)(3,l)  (1,3)  (3,  x)  (!,«)(«,  l)  (1,8)  (8,1)  (l,l)  (1,1)  (L,l) 

in  which 


1 


(43) 

gfo)  =  [T  ra  r  r  f(E) ] ,  =  (ax  afl  a3  a4 


BW  =  (1),  6«  =  (^). 


Before  proceeding,  we  shall  address  one  possible  problem  with  Geoceiver  tracking  that  ha* 
not  so  far  been  considered.  This  is  concerned  with  what  to  do  in  the  case  of  one  or  more 
temporary  tracking  dropouts  during  a  (.ass.  When  phase  lock  is  lost,  the  doppler  counter 
is  immediately  set  to  zero  and  remains  set  to  zero  until  the  first  one  minute  mark  after 
phase  lock  is  restored.  Thus  a  zero  cycle  count  for  a  given  counting  interval  is  a  positive 
indication  of  a  dropout  during  that  interval.  When  counting  is  resumed  following  a 
dropout,a  new  constant  of  integration  r0,  or  equivalently  a  new  error  coefficient  a0,  must 
be  established.  On  the  other  hand,  it  is  clear  from  their  physical  interpretation  that  the 
remaining  error  coefficients  (05,03,03,04,05)  will  not  be  altered  by  a  dropout.  Hence, 
only  a0  requires  reinitialization  following  a  dropout.  For  each  pass  provisions  are  made 
in  SAGA  to  accommodate  up  to  a  maximum  of  three  dropouts  in  tracking  from  each 
participating  station.  This  is  accomplished  by  writing  the  observational  equation 
generated  by  the  jth  point  as: 

(44)  Aj  Vj  +  Bjfi=  <, 

in  which  • 
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(45) 

Bt  = 

-c 

rM  o(3*)  B(a> 

°P i  °P3  °PJ 

SlpB»j 

t  fti4) 
Wj  BS3 

«».  B(p! 

t  dW 

V*p  BP1 

(1,13+4) 

(l/2>) 

(i,a)  (i,e)  (l, b) 

(l/l) 

(l/l) 

(I'X) 

(l/l) 

(46) 

5T 

(  5<1>T 

& 

6<*>V 

(1,13+t) 

(l,3) 

(ir a)  fit  ® )  (i/*) 

(l,l) 

(i/i) 

(i/i) 

(l/l) 

in  which  the  subscript  p  connotes  the  pth  tracking  interval  (max  p=4)  and 

(47) 


«t,  =  l  ^  i “P 
4ip  =0  'f  i^P  • 


The  dimension  6  indicates  the  number  of  error  coefficients  exercised  by  the  station  on  the 
pass  (this  can  range  from  a  minimum  of  six  to  a  maximum  of  nine  depending  on  the  number 
of  tracking  dropouts). 


3.4  Normal  Equations 

The  development  of  the  normal  equations  gv..nerated  by  Geoceiver  observations 
follows  precisely  the  process  outlined  in  Section  2.7  for  optical  observations.  With 
suitable  (and  perfectly  obvious)  reinterpretation/equations  (44)  through  (49)  hold  also 
for  Geoceiver  observations.  In  particular,  the  partial  set  of  normal  equations  generated 
by  a  given  station  for  a  given  pass  (eq.  (49)  of  Sec.  2.7)  applies  equally  to  Geoceiver 
observations.  So  also  does  the  merged  system  of  equations  from  all  stations  for  a  given 
pass  (eq.  (50),  Sec.  2.7),  with  one  proviso,  namely,  that  each  Csoceiver  station  be 
assigned  a  distinct  station  number  even  when  it  is  colocated  with  an  optical  station. 
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This  step  automatically  accounts  for  the  fact  that  the  vector  6^  of  error  coefficients 
for  a  Geoceiver  is  different  from  the  vector  fift  of  error  coefficients  for  a  camera.  Should 
a  Geoceiver  be  colocated  with  a  camera,  appropriate  interstation  constraints  can  be 
exercised  to  insure  that  both  stations  receive  a  common  adjustment.  With  this  under¬ 
standing  then,  we  can  regard  the  general  normal  equations  (52)  as  applying  equally 
well  to  (a)  a  network  of  optical  tracking,  (b)  a  network  of  Geoceivers  (or  other 
ranging  systems),  (c)  a  combined  network  of  optical  trackers  and  Geoceivers. 

3 . 5  A  Priori  Constraints 

Of  the  six  coefficients  ir.  the  Geoceiver  error  model,  the  first  two  a0  and  ax  have 
large  effects  and  are  unlikely  to  be  subject  to  worthwhile  a  priori  constraints.  The  remaining 
four,  on  the  other  hand,  have  relatively  small  effects  and  are  subject  to  rather  tight  a  priori 
constraints.  When  the  Geoceiver  is  functioning  properly  and  when  the  frequency  bias  6f0 
of  the  satellite  oscillator  is  reasonably  well  monitored,  the  a  priori  constraints  indicated 
in  the  table  below  can  serve  as  rough  guidelines. 


Coefficient 

Physical  Significance 

Typical  A  Priori  Constraint 
(One  Sigma) 

«o 

(i+f), 

'  *00  ' 

10s  to  107m 

ai 

10,\o  to  TO3  X0  m/sec 

°3 

|X0(f0-fi) 

Xofox10~n  m/sec  s 

a3 

-6f0/fco 

-e  ,  •ft 

10  to  10 

°4 

fir  (timing  bias) 

50fi  s 

ap 

fin;  (refraction) 

0.2  m 

In  our  view,  the  major  shortcoming  of  Geoceiver  is  its  low  sampling  rate  of  only 
one  readout  per  minute.  This  means  that  many  passes  will  yield  under  ten  observations 


and  even  the  longest  passes  are  unlikely  to  yield  more  than  twenty  observations.  The 
justification  of  the  low  sampling  rate  probably  stems  from  the  use  (in  one  mode  of  operation) 
of  Transit  timing  signals  to  trigger  readout  of  cycle  count.  Also,  a  probable  factor  is  the 
approach  to  data  reduction  envisioned  by  the  designers  of  Geoceiver.  As  outlined  in 
Stansell,  et.  al.  (1965),  this  approach  involves  using  an  observational  equation  of  the 
form: 

(48)  r3  -  =  X0  ANj  +  (f0  -  fo)(r  3  -  T 3-i)  +  neglected  higher  order  terms, 

which  is  relatively  weak,  geometrically,  for  small  time  intervals.  This  formulation  has  the 
advantage  of  completely  eliminating  the  error  coefficient  for  zero  set  a0  and  of  generally 
suppressing  to  insignificance  the  effects  of  a2/a3,a4,a5  (over  intervals  as  short  as  one 
minute).  On  the  other  hand,  it  fails  to  exploit  one  of  the  strongest  characteristics  of  the 
Geoceiver,  namely,  the  nondestructive  readout  of  cycle  counts.  It  is  clear  that  (48)  is 
equally  valid  for  destructive  or  nondestructive  readout  of  cycle  counts,  whereas  our 
formulation  (eq.  (131))  is  designed  specifically  to  exploit  the  nondestructive  character  of 
the  readout.  If  readout  is  indeed  nondestructive,  our  approach  should  result  in  considerably 
more  effective  utilization  of  Geoceiver  observations.  On  the  other  hand,  our  approach 
could  benefit  considerably  from  an  increased  sampling  rate .  An  internally  triggered  readout 
every  ten  or  twenty  seconds  would  probably  be  close  to  ideal .  Such  serial  correlation  as 
might  thereby  be  introduced  could  be  properly  taken  into  account  by  the  process  of  auto¬ 
regressive  feedback  developed  in  Section  5. 

With  the  readout  rate  as  low  as  it  is,  one  cannot  generally  expect  to  obtain  from 
the  adjustment  any  significant  degree  of  improvement  (over  a  priori  values)  in  the  adjusted 
error  coefficients  a3,a4,a5.  Here,  what  one  mainly  accomplishes  is  the  realization  of  a 
more  valid  error  propagation  by  considering  the  effects  of  these  sources  of  error.  The 
a  priori  value  of  a3  is  subject  to  a  degree  of  improvement  ranging  from  slight  for  short 
passes  (under  5  minutes)  to  pronounced  for  moderately  long  passes  (over  15  minutes).  On 
long  passes  observed  by  several  starions  (five  or  more),  one  can  expect  to  determine  values 
of  ax  and  as  to  accuracies  (sigma)  of  better  than  five  meters  and  five  millimeters  per 
second,  respectively. 
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3.6  Reduction  of  Other  Ranging  Systems 

Although  emphasis  in  this  section  has  been  on  the  Geoceiver,  our  treatment  of 
the  Geoceiver  as  a  ranging  system  makes  it  possible  to  employ  SAGA  for  the  reduction 
of  ranging  systems  in  general.  The  primary  requirement  for  processing  other  ranging 
observations  through  SAGA  is  that  they  be  presented  to  the  program  in  the  Geos  data 
format.  With  active  ranging  systems,  the  coefficient  a0  would  ordinarily  be  subject  to 
very  tight  a  priori  constraints  (e.g.,  a  few  meters  for  lasers),  and  the  coefficients  alf 
a2  would  be  inapplicable  and  hence  would  be  suppressed.  An  option  is  provided 
within  SAGA  for  the  application  of  corrections  for  tropospheric  refraction  In  the  event 
that  such  corrections  are  lacking. 
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4.0  THEORY  AND  EXECUTION  PARTITIONED  REGRESSION 


4. 1  First-  Older  Partitioned  Regression 

SAGA  is  designed  to  process  an  arbitariiy  large  number  of  passes,  each  observed 
by  a  subset  of  as  many  as  fifteen  stations  each  of  which  can  conceivably  (with  re¬ 
initialization  of  error  coefficients)  introduce  as  many  as  seventeen  error  coefficients 
peculiar  to  the  pass.  Accordingly,  SAGA  generates  a  patterned  system  of  normal 
equations  that  grows  in  dimensions  both  with  the  number  of  passes  processed  and  with 
the  number  of  error  parameters  exercised  on  each  pass.  In  practical  application:;,  the 
normal  equations  can  grow  to  embrace  thousands  of  unknowns  and  thus  be  unamenable 
to  solution  by  conventional  reductions  that  fail  to  exploit  the  patterned  characteristics 
of  the  system.  In  this  section  we  shall  develop  practical  agorithms  for  the  solution  of 
normal  equations  related  to  those  generated  by  SAGA. 

The  normal  equations  for  a  general  regression  analysis  may  be  written  as 

(i)  N6  “  c 

in  which  6  denotes  the  parametric  vector  to  be  estimated,  N  is  the  coefficient  matrix, 
and  c  is  the  constant  column.  Let  6  be  arbitrarily  be  partitioned  into  two  vectors 
6,  6  so  that 


(2)  6  = 


Then  the  original  normal  equations  partitioned  to  be  conformable  with  the  partitioning  of 
6  can  be  written  as 

(3)  [n  nIIYI  I’d 

NT  N  6  c  . 

To  this  point  the  normal  equations  are  of  perfectly  general  form.  We  now  abandon  full 
generality  by  assuming  that  the  parametric  vector  6  is  common  to  most  if  not  all,  of  the 
original,  observational  equations,  whereas  6  is  composed  of  (possibly)  a  large  number  of 
subvectors  *6^  ,  <S3, . .  .,'6B,  each  of  which  appears  only  in  a  subgroup  of  observational 
equations  to  the  exclusion  of  other  subvectors  of  C.  Under  these  circumstances,  the 
normal  equations  (3)  assume  the  specialized  form 


N  N,  ...  N,  6  c 

___! - - — - 

Rf|f^  o  ...  o  £  <4 

i 

(4)  N£  i  0  N,  . . .  0  6a  =  ca 

*  I  •  •  •  • 

:  j :  :  _ o  :  : 

Rt I  0  0  ...  N,  fi  C,  . 

Here,  the  N  portion  of  the  matrix  assumes  a  block  diagonal  form  by  virtue  of  the  assumption 
that  for  all  i  /  j,  and  flj  do  not  appear  in  common  observational  equations. 

To  derive  this  result  we  may  proceed  as  follows.  The  linearized  observational 
equations  may  be  oi  Jered  into  groups  corresponding  to  the  various  parametric  subvectors 
6,.  Thus  we  write 

v,  +  ^  fi  +  \ 

(5)  va  +  Ba  6  +  B,  5a  = 

vn  +  &n  6  +  Bn  6n  " 

•  •  • 

Here  the  v's  ore  residual  vectors,  the  Bj's  and  Bj's  are  the  coefficient  matrices  of  the 
*  *  * 

parametric  vectors  6  and  6*#  respectively,  and  the  are  the  discrepancies  between 
actual  observations  and  their  computed  values.  We  note  that,  in  accord  with  the  above 
discussion,  the  6  vector  is  distinguished  by  being  common  to  all  groups  of  the  observational 
equations,  whereas  each  *5,  vector  appears  in  one  and  only  one  group  of  observational 
equations.  The  above  system  may  be  written 


vx  Bj  0  . . .  0 

6x 

(6)  va  +  Ba  0  Ba  ...  0 

•  ■  •  •  •  V 

«  •  •  *  <  Oq 


vj  [_Ba  0  0  ...Bj  .* 


1 


or  more  compact!  •  as, 

(7)  v  +  3  5  -  c  . 


If  A  denotes  the  covariance  matrix  of  the  vector  of  observations,  the  minimum  variance 
adjustment  is  obtained  determining,  of  vectors  v  and  6  that  satisfy  (7),  that  particular 
pair  that  leads  to  the  minimization  of  the  quadratic  form 

(8)  s  =  vTA_1v. 


This  process  leads  to  the  normal  equations 

(9)  Nfi  =  c. 


where 

(9a)  N  =  BT  A"1  B  / 

(9b)  c  =  BT  A"1  €  • 

We  shall  now  assume  that  A  is  of  the  form 

(10)  A  =diag  (Ax  Aa  •  •  •  A.) 

where  each  Ai  consists  to  the  covariance  matrix  of  the  observational  vector  corresponding 
to  the  residual  vector  v4 .  Then  inasmuch  as  B  and  e  are  of  the  form* 


(ii)  b  =  :  :  : 


6^  Bj  0  . . .  0 

Ba  0  Ba  ...  0 


B  0  0  . . .  B. 

u  0 


r  €  “  * 


It  follows  by  direct  evaluation  of  (9a  )  and  (9b  )  that  the  normal  equations  (9)  are  of  the 
form  (4)  in  which 

(12)  N-=  Nj_  +  N3  +  . . .  +  N#,  c  *  Cj  +  c„  +  . . .  +  cB 
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wherein 


I 

r ' 


!i 


it 

[; 

t 


-1 


03)  Nt  =  B'l  At  B,, 


A?  €, 


and  in  which 

(u)  Nt = bt  A;x 

05)  Nj  =  bt  a;1  &V  c.  =bt  a;1  c,  . 


It  is  to  be  noted  that  this  result  applies  to  any  set  of  observational  equations  of  the  form  (5) 
for  which  the  covariance  matrix  of  the  observational  vector  is  of  the  form  (10).  Systems  of 
this  nature  arise  in  a  broad  range  of  applications.  In  particular,  such  systems  are  generated 
in  trajectory  and  orbital  analysis  when  coefficients  of  tracking  error  models  are  to  be 
recovered  simultaneously  with  trajectory  p  rameters.  Before  taking  up  specific  examples, 
however,  we  shall  continue  in  cn  abstract  vein  to  study  the  properties  of  systems  of  nornial 
equarions  having  the  particular  patterned  structure  of  (4). 

Returning  then,  to  equation  (4),  we  first  remark  that  when  the  number  of  elements 

e  •  « 

in  the  fi  vector  is  small  compared  with  the  number  in  the  6  vector,  the  block  diagonality 
of  the  N  portion  of  the  matrix  is  but  of  minor  advantage  to  the  solution  of  the  system. 
However,  in  many  problems  the  order  of  5  is  much  greater  than  that  of  6,  while  at  the 

*  e  e 

same  time  the  order  of  each  6,  is  either  smaller  than  that  of  6  or  else  is  of  roughly  comparable 
magnitude.  Moreover,  the  number  n  of  5j  submatrices  may  be  very  large  and  without  any 
particular  preset  limit.  Here,  the  normal  equations  can  assume  enormous  dimensions  and 
the  possibility  of  a  practical  solution  hinges  entirely  on  the  block  diagonality  of  Kl. 

We  shall  call  a  system  of  normal  equations  of  the  form  (4),  a  first  order  partitioned 
system  and  shall  express  the  coefficient  matrix  schematically  as 


06) 


*54" 


l 


l 

\ 


1 


* 


! 


The  sense  of  this  dicgram  is  the  following:  the  horizontal  an d  vertical  arrows! represent 
solid  rows  and  columns,  respectively,  of  elements  that  are  not  necessarily  zeijo;  the 
diagonal  arrow  represents  a  block  diagonal  array  of  submafrires;  the  blank  regions 
between  the  arrows  represent  the  portions  of  the  matrix  thot  are  completely  filled  with 
zero  elements.  The  arrowheads  are  indicative  of  the  fact  that  the  size  of  the  matrix  can 
grow  without  a  preset  limit.  We  shall  find  this  diagrammatic  representation  of  a  first 
order  partitioned  system  to  be  a  convenient  point  of  departure  in  our  later  discussions  of 
higher  order  partitioned  systems. 


Aside  from  noting  that  the  block  diagonaiity  of  a  first  order  partitioned  system 
provides  the  key  to  its  practical  solution,  we  have  vet  to  develop  the  details  of  the 
solution.  We  shall  now  remedy  this  deficiency.  To  begin,  let  us  formolly  define  the 
inverse  of  N  to  be 


(17) 


M  M 

MT  M 


t  _  *'-•«  V  tmm  •• 

where  M,  M,  M  are  each  of  the  same  order,  respectively,  as  their  counterparts  N,  N,  N 
in  N.  Inasmuch  as  NM  =  I,  the  unit  matrix,  by  virtue  of  the  definition  of  M,  we  may 
write  the  equivalent  relation  in  partitioned  form  to  get 


1 - 

IZ 

z 

1 _ 

M  M 

r,  a] 

|_Nr  N_ 

MT  M 
■*  - 

- 1 

o 

_ 1 

Performing  the  multiplication  on  the  left  hand  side  of  the  equation  and  equating 
the  resulting  elements  with  the  corresponding  elements  on  the  right,  we  arrive  at  the 
foliowing  four  matrix  equations: 


(19a)  NM  +  NM*  =  I 
(19b)  NM  +NM  =0 
(19c)  NTM  +  NMT  =0 
(19d)  N'M  +  NM  =1  . 
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Since  we  have  already  exploited  the  fact  that  the  inverse  of  a  symmetric  matrix 
is  itself  symmetric,  the  above  system  of  four  matrix  equations  involves  but  three 
unknown  matrices;  M,  M,  M.  Various  mutually  consistent  solutions  are  possible 
depending  on  which  set  of  three  one  elects  to  solve.  In  view  of  the  fact  that  we  are 
ultimately  going  to  regard  N  os  large,  but  block  diagonal,  it  turns  out  appropriate  to 
select  equations  (19a),  (19c),  and  (19d)  for  the  solution.  Accordingly,  first  solving 
(19c)  for  Mt  in  terms  of  M,  we  get 

(20)  MT  =  -N_1  NT  M 

and  upon  substituting  this  into  (19a)  and  (19c),  we  have 
(21a)  NM  -  NN'1  NT /v\  =  I, 

(21b)  ~NTMNN‘1  +  NM  - 1. 

Thu  solution  of  (21a)  for  M  and  that  of  (21b)  for  M  gives 
(22a)  M  »  (N  -  RN"lRT)“X 
(22b)  M  =  N-t+  N-1  NT  M  N  Kf 1  . 

Once  M  has  been  obtained  from  (22a),  the  result  can  be  substituted  in  (22b)  to  obtain 
M  and  in  (20)  to  obtain  MT  .  Thereupon,  the  formal  solution  of  the  original  normal 
equations  (2)  may  be  written 


(23)  6  = 


N  N 


RT  N 


M  M 


MT  M 


(24a)  6  =  M  d  +  Me, 

(24b)  *5  -•  MT  c+  Me. 

Alternatively,  by  virtue  of  the  expression  for  M  given  by  (20),  6  can  be  written  as 
(25)  6  =  M(c  -  R  N  d) . 
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Similarly,  by  virtue  of  (20)  and  (22b)  the  expression  (24b)  for  '6  becomes 

(26)  a  =  -N"1  NT  M  6  +  (N^c  +  n”1  Nt  M  N  IsJ'Vc 
The  right  hand  side  of  which  can  be  factored  to  yield 

(27)  6  =  N"1  c  -  N”1  Rt  M(c  -  R  N"1  c)  . 

Equation  (25)  permits  this  to  be  written  in  the  more  compact  form 

(28)  6  =  N_1  c  -  N"lNT  fi. 


Equations  (22a),  (25)  and  (28)  constitute  the  major  preliminary  result*  that  we  seek. 
They  are  preliminary  because  we  have  yet  <o  exploit  the  block  diagonality  of  Kl.  ThU 
immediately  permits  us  to  write 


(29)  N_1  = 


••  —1 

rV  o 


...  o 


Na  ...  0 


...  0 


_  ■ 
0  ...  N 


-x 


Since  N  is  of  the  form 

(30)  N  =  (Nj  Na  ...  NJ, 

the  expression  Q  =  N  1  NT  becomes 


(31)  Q  = 


N"1  Ni“ 

V 

K  n; 

II 

q8 

_n;£  % 

—  «■ 

(32) 


and  the  expression  R  =  N  N  1  RT  =  RQ  becomes 

R  ■  Nx  N”1  NJ  +  N8  +  -  •  +  NB  N”1  RTa 
=  NiQa  +N8Qa +  ...  +RaQa 
—  Rj  +  Ra  +  . . .  +  R 
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"'  •|]>iilb  Ifciimwnffc1rt,il 


rhr  r  atOMtatmam 


-  .rfi  fcii  1  rMsilMhe  r  inrf  i  eiMfis  t  iMu  sfiaefl  .ttArviTi  -fiikfa 


At  this  point  we  recall  from  equation  (12)  that  the  matrices  N  and  c  can  be  expressed  as 
(33a)  N  =  Nj  +  Na  +  . . .  +  N# 

(33b/  c  =  Cj  +  ca  +  . . .  +  cn 

in  which  N,,  c,  are  generated  by  the  particular  subset  of  observational  equations 
containing  5,  (and  hence  by  the  same  subset  that  generates  Nt,  N„  and  c,).  It  follows 
then,  that  M  can  be  written  as 

(34)  M  =  Nj  +  N#  +  .  * .  +  N„  -  (\  +  Rs  +  . . .  +  R„) 


(35)  M  =  Sj  +  Sa  +  . . .  +  S, 
where 

(3b)  S,  =  -  Rt . 

Similarity,  the  expression  6  -  NN  1’c  in  (25)  assumes  the  form 

(37)  c  -  N  N_1c  =  («j  +  cs  +  . . .  +  ca)  -  (Q^c*  +  Q8c3  +  . . .  +  Qa  c.) 


=  +  . . .  +  r8 


where 


(38)  Cj  c,  “  Q,c,. 

Proceeding  further  to  examine  the  consequences  of  the  block  diagonality  of  the  N  portion 
of  the  normal  equations,  we  find  that  the  expression  (28)  for  "fi  can  be  written  as 

fij  Njl  0  ...  0  cx  Qi 

(39)  6=  6a  =  0  ...  0  c8  -  Qs  6 


o  o  ...Mr  c 


i*wnew.iwn 

'■  - 


From  this  if  follows  thaf 

(40)  «,  =  N"1  c,  -  Q  t  5, 

which  shows  that  each  vector  can  be  independently  determined  once  the  value  of  £ 
has  been  obtained. 


Inasmuch  as  the  inverse  of  the  coefficient  matrix  of  the  normal  equations  provides 
the  covariance  matrix  of  the  adjusted  parameters,  it  follows  immediately  from  the 
foregoing  development  that 

(41)  var  (6)=/0\. 

Since  M  has  to  be  computed  in  order  to  obtain  the  solution  for  fi,  the  covariance  matrix  of 
6  falls  out  as  a  direct  by-product  of  the  solution.  A  corresponding  result  does  not  hold  for 
the  vector  5,  for  although  the  covariunce  matrix  of  *6  is  given  by  M,  we  found  it  possible 
to  by  .  oss  the  evaluation  of  M  in  the  computation  of  6  by  equation  (40).  If,  however,  we 
return  equation  (22b)  and  trace  through  the  consequences  of  the  block  diagonality  of  N, 
we  shall  find  that  the  diagonal  submatrix  of  M  defining  the  covariance  matrix  of  fl,  is 
given  by 

(42)  var  (6t)  =  N«  +  Q  s  M  QT, . 

Furthermore,  the  covariance  matrix  of  an  arbitrary  pair  of  vectors  6,,  is  given  by 

(43)  cov  (61,54)=Q,MQt1,  (i/j). 

This  completes  the  derivational  development  of  the  theory  underlying  first  order 
patterned  regression.  Because  the  essential  computational  flow  is  perhaps  obscured  by 
derivational  detail,  it  is  appropriate  that  we  conclude  this  section  by  extracting  the 
essentials  of  the  reducNon.  We  shall  take  as  our  starting  point  the  following  set  of  four 
matrices  concerned  with  the  i  th  group  of  observational  equations:  B{,  B4,  At. 

Starting  with  i  =  1,  one  first  computes  in  terms  of  these  matrices,  the  following  five  matrices: 


I 


■»»«« jaMt*  -■‘traiitivi  .iriiftiiH 


1 


(44) 

N. 

=  b; 

Al'Bt 

(45) 

N| 

a;1  b\ 

(46) 

Nt 

-  K1 

Al1  Bt 

(48) 

A^  £, 

(49) 

«  • 
ci 

=b; 

-X 

A,  c, 

In  terms  of 

these 

the  follow 

(50) 

Q. 

=  n; 

lNf 

(51) 

R« 

=  Rt 

Q. 

(52) 

st 

=  N, 

(53) 

• 

=  ci 

-QT.*e*i 

As  St  and  ct  are  computed,  they  are  added  to  the  sum  of  their  predecessors  and  only  the 
cumulative  sum  is  retained.  After  all  groups  of  equations  have  been  processed  sequentially 
in  this  manner,  one  arrives  at  the  final  values. 

(54)  S  =  Sj  +  S0  +  . . .  +  SB , 

(55)  c=^  +?8  +  ...  +?„ 

The  solution  for  6  is  then  obtained  immediately  from 

(56)  6  =mT  j 

•  “1  *  •• 
in  v/nich  M  =  S  is  also  the  covarianGe  matrix  of  fi.  Thereupon  the  solution  for  each  fit, 

in  turn,  is  computed  from 

(57)  6,  -  Nj1  ct  -  Q  t  fi  , 

parallel  with  which  the  covariance  matrix  is  computed  from 

,  _  _ 


(58)  .var  (i,)  =  Klj  +QtMQ][. 


Final  I/,  the  observational  residuals  are  obtained  by  substituting  the  values  determined 
for  6  and  5,  into  the  original  observational  equations,  which  gives 

(60)  vt  -  -  Bj  6  -  8j  6,  • 

The  salient  properties  of  the  above  solution  of  the  normal  equations 
characterizing  first  order  patterned  regression  are 

(1)  The  overall  computational  effort  generally  tends  to  increase  oniy  linearly 
with  the  number  of  parameters  being  recovered,  rather  than  with  the  cube  of 
the  number  as  in  adjustments  where  patterning  is  not  or  cannot  be  exploited. 

This  renders  practical  the  solution  of  important  problems  involving  the 
simultaneous  determination  of  thousands  or  tens  of  thousands  of  unknowns. 

(2)  The  process  of  solving  the  normal  equations  proceeds  apace  with  the 
formation  of  the  normal  equations  Thus  by  the  time  the  last  group  of 
observational  equations  ha:  b.en  processed,  the  great  bulk  of  the  computation 
required  for  the  solution  has  also  been  completed. 

(3)  Core  storage  requirements  are  mininjal  in  a  computer  program,  since 
the  processing  for  each  of  the  basic  groups  of  observational  equations  can 
proceed  independently  of  that  of  the  other  groups. 

Before  taking  up  some  of  the  applications  of  first  order  partitioned  regression,  we  would 
note  that  the  algorithm  (equations  (44)  to  (53))  leading  to  the  computation  of  5t,  <?,  can 
be  put  into  the  alternative,  more  compact  form 


.  his  representation  facilitates  comparison  with  the  algorithm  to  be  developed  for  second 
order  partitioned  systems.  The  application  of  first  order  partitioned  regression  would  have 
been  sufficient  for  SAGA  were  it  not  for  the  feet  that  the  program'*  general  capabilitie* 

•  e 

with  respect  to  error  modeling  can  lead  to  very  large  Nt  matrices. 
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4.2  Applications  of  First  Order  Partitioned  Regression 

The  original  application  of  the  above  development  was  made  in  Brown  (1958a), 
where  it  provided  the  basis  for  the  adjustment  of  a  genera!  photogrammetric  net.  Here,  the 
6  vector  corresponded  to  the  unknown  elements  of  orientation  of  the  cameras  and  each  6, 
vector  corresponded  to  the  unknown  X,Y,  Z  coordinates  of  a  photographed  point.  One 
specialized  application  of  this  theory  was  to  satellite  geodesy  (Brown,  1958b),  wherein 
the  6  vector  consisted  of  the  unknown  coordinates  of  a  set  of  cameras  and  the  <5,  vectors 
consisted  of  the  unknown  X,Y,Z  coordinates  of  recorded  flashes.  This  theory  was 
successfully  applied  to  the  determination  of  the  precise  location  of  the  location  of  Bermuda 
relative  to  the  North  American  continent  (Brawn,  1959)  and  later  to  a  detailed  study  of 
the  feasibility  of  employing  similar  techniques  on  a  larger  scale  to  produce  an  ultra-precise 
survey  of  key  tracking  stations  along  the  ’  . Icntic  Missile  Range  (Brown,  1960). 

Another  class  of  applications  of  first  order  partitioned  regression  emerged  in  1959, 
when  it  was  recognized  that  complex  problems  of  trajectory  analysis  could  be  formulated 
in  a  manner  leading  ta  normal  equations  of  essentially  the  same  form  as  those  successfully 
handled  in  the  photogrammetric  application.  Here,  the  fi  vector  consisted  of  unknown 
coefficients  of  tracking  error  models,  whereas  the  *5,  consisted  of  unknown  X,Y,Z 
coordinates  of  trajectory  points.  This  application  (Brown,  et.al.1961)  became  known  as 
EMBET  (Error  Model  Best  Estimate  of  Trcjectory). 

In  due  course  it  became  recognized  through  application  of  EMBET  to  conventional 
tracking  systems  that  certain  systematic  errors  (e.g.,  zero  set  errors)  could  often  be 
recovered  more  accurately  through  data  reduction  than  they  could  be  established  by  means 
of  hardware.  In  1960  Brown  suggested  that  self-calibration  by  means  of  EMBET  should  be 
exploited  as  a  guiding  principle  in  the  very  design  of  new  tracking  systems.  This 
philosophy  was  adopted  in  the  design  of  the  GLOTRAC  system.  Here,  unknown 
reference  frequencies  of  remote  tracking  stations  were  recovered  in  a  specific  EMBET 
reduction  called  GLAD  (GLotrac  ADjustment). 


-62- 


In  Brown  (1964b)  and  in  Brcwn  et.  al ,  (1963),  (1964),  the  feasibility  of  self- 
calibration  of  tracking  systems  by  means  of  observation*  of  satellites  was  extensively 
studied  and  simulated.  This  led,  in  particular,  to  an  important  version  of  EMBET 
called  NEO-EMBET  (N  Epoch  Orbital  -  Error  Model  Best  Estimate  of  Trajectory) 
NEO-EMBET  renders  practical  the  simultaneous  reduction  of  observations  of  an  unlimited 
number  of  satellite  arcs  with  all  arcs  being  interrelated  by  certain  common  parameters 
(such  as  coordinates  of  tracking  stations  or  stable  coefficients  of  error  models)  but  with 
each  arc  requiring  the  recovery  of  a  fresh  set  of  orbital  elements  and,  in  addition,  with 
each  pass  observed  by  a  given  tracker  requiring  recovery  of  reinitialized  error  coefficients# 

In  a  study  performed  for  NASA  (Brown,  Stephenson,  Hartwell,  1965)  it  was 
recommended  that  a  special  NEO-EMBET  reduction  be  implemented  in  support  of  the 
GEOS  I  Short  Arc  Experiment.  This  recommendation  was  adopted  by  NASA  and  led  to 
the  development  of  an  unusually  powerful  ^  -id  flexible  program  called  GDAP  (GEOS 
Adjust,  .  ,,t  Program)  which  is  c!i  ussud  by  Brown  (1967a)  and  is  described  in  detail  by 
Lyri.-i  , '  967).  A  significant  application  of  GDAP  to  the  establishment  of  a  much 
improved  survey  of  the  GEOS  short  arc  tracking  network  is  reported  by  Brown  (1968). 

Under  a  recently  completed  contract  with  NASA,  GDAP  was  employed  as  the 
starting  point  for  the  development  of  a  still  more  advanced  program  called  NAP 
(Network  Analysis  Program).  Whereas,  GDAP  can  accomodate  an  unlimited  number  of 
interrelated  short  arcs,  NAP  is  able  to  accomodate  an  unlimited  number  of  interrelated 
long  arcs.  However,  the  long  arc  application  introduces  what  we  shall  presently  take 
up  as  second  order  partitioned  regression  As  has  already  been  indicated,  the  purpose 
of  the  first  part  of  this  section  is  to  provide  the  background  needed  for  an  understanding 
of  second  order  partitioned  regression  Before  taking  up  this  topic,  we  would  briefly 
cite  rive  more  instances  where  first  order  partitioned  regression  has  been  successfully 
applied. 

In  Brown  (1964),  the  technique  was  applied  to  the  reduction  of  stellar  plates 
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recorded  by  metric  cameras  of  long  focal  length.  Here,  random  errors  in  cataloged 
stellar  positions  become  significant.  Accordingly,  the  reduction  was  formulated  so 
that  parameters  peculiar  to  the  camera  (focal  length,  distortion,  orientation,  etc.) 
were  incorporated  into  the  6  vector,  whereas  the  corrections  to  the  cataloged  positions 
of  the  I  th  star,  were  accounted  for  by  6 1 . 

In  Brown  (1966)  EMBET  and  NEO-EMBET  reductions  were  developed  and 
successfully  applied  to  the  reduction  of  Geodetic  SECOR  observations. 

In  Brown  and  Trotter  (1967),  first  order  partitioned  regression  provided  the  basis 
for  the  development  of  the  Method  of  Continuous  Traces,  a  technique  for  exploiting 
measurements  of  uninterrupted  photographic  traces  of  sun-illumined,  passive  satellites 
to  establish  precise  geodetic  positions.  Tire  technique  does  away  with  the  need  for 
synchronized  chopping  shutters,  thu.  eadmg  to  less  expensive,  much  simplified,  and 
more  reliable  data  acquisition. 

In  Brown  (1967b)  the  development  and  successful  application  of  a  plate  measuring 
comparator  is  described.  The  design  of  the  comparator  is  based  directly  on  principles  of 
self-calibration  as  made  practical  through  first  order  partitioned  regression. 

In  Brown  (1969)  the  calibration  of  metric  cameras  was  approached  from  the  stand¬ 
point  that  parameters  of  the  inner  cone  (radial  and  deconturing  distortion,  principal 
distance,  and  coordinates  of  principal  point)  are  invariant  over  an  indefinitely  large 
number  of  exposures,  whereas  elements  of  exterior  orientation  may  change  from  exposure 
to  exposure.  This  formulation  of  the  problem  of  camera  calibration  led  to  the  formation 
of  a  first  order  partitioned  system  of  normal  equations. 

From  the  foregoing  discussions  it  is  amply  clear  that  the  theory  of  first  order 
partitioned  regression  has  been  put  to  extensive  use  during  the  course  of  the  last  decade. 
Undoubtedly  many  more  applications  will  emerge  in  the  next  few  years.  There  are, 
however,  as  in  SAGA,  significant  problems  leading  to  highly  patterned  normal  equations 
that  are  not  amenable  to  solution  by  the  first  order  scheme.  Many  such  problems  fall  into 
the  province  of  second  order  partitioned  regression  which  is  the  topic  of  the  remainder  of 
this  section.  We  expect  that  during  the  next  decade  a  host  of  applications  of  second 
order  partitioned  regression  will  emerge  in  various  disciplines. 
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4. 3  Second  Order  Partitioned  Regression 

As  in  our  treatment  of  first  order  partitioned  regression,  we  shall  proceed  first  in 
an  abstract  vein  in  developing  the  theory  of  second  order  partitioned  regression.  We 
begin  by  postulating  a  patterned  system  of  normal  equations  of  the  form 

A  rv;  rw  fW  ""1  P/y  "1  f** a 

U  Ut  Ua  . .  .  UB  6  j  C 

UJ  £4  0  ...  0  6i  ct 

(61)  TJJ  0  CJ8  ...  0  0a  =  Cg 


\?i  0  0  ■  •  * a J  L^-J  LcJ 

which,  ir  will  be  noted,  has  the  .ame  structure  as  the  first  order  system  defined  by 
equc.. '.on  (4).  Also,  as  in  (4),  the  index  n  can  be  indefinitely  large.  What  distinguishes 

rw  m 

(61)  as  a  second  order  system  is  its  finer  structure.  The  submatriees  U1#  Utr  5t  end  ct  are 

of  the  form 


(62a)  Ut  =  [ut  Usl  Ut3  ... 


N{  Nu  Nja  . . .  N1Bj 

NTu  Nji  0  ...  0 


(62b)  *Uj  -  Nja  ®  Nta  ...  0  ,  6j  -  6ia  ,  ct  =  c\s 


nt1Bjo 


0  ...  N, 


in  which  the  index  n4  can  be  indefinitely  large.  It  is  particularly  to  be  noted  that  each 
U  has  the  structure  of  a  first  order  system  (and  hence  can  grow  without  preset  limit).  This 
coupled  with  the  fact  that  the  number  of  U,  is  also  without  preset  limit  is  the  essence  of 
a  second  order  patterned  system.  Such  systems  may  be  represented  diagrammatical ly  os  in 


Mr  rrnnrmnrri  r  w  ■ 


Figure  4,  Here,  the  long  horizontal  and  vertical  arrows  represent  solid  rows  and 
columns  of  nonzero  elements  and  the  blocks  arrayed  along  the  diagonal  individually 

consist  of  first  order  patterned  systems.  Thus  a  second  order  system  may  be  described 

{ 

as  a  first  order  system  of  first  order  systems.  Similarly,  a  third  order  system  may  be 
described  as  a  first  order  system  of  second  order  systems.  Although  we  shall  not  be 
directly  concerned  with  the  reduction  of  third  and  higher  order  systems,  we  have 
provided  a  diagrammatic  representation  of  a  third  order  system  in  Figure  5.  While  we 
have  yet  to  find  a  solid,  practical  application  for  a  third  order  scheme,  there  do 
exist,  as  we  have  already  suggested,  many  significant  applications  of  second  order 
schemes. 


Let  us  proceed  formally  with  (61)  as  if  it  were  a  first  order  system.  We  could 

then  follow  the  procedure  of  Section  1  to  ootain  a  solution  first  for  6  and  then  for  each 

of  the  *6  ,  in  turn.  The  practical  difficulty  with  this  approach  is  that  the  *U,>  unlike 
•  • 

their  counterparts  N,  in  a  first  order  scheme,  may  grow  without  set  limit.  While  it  is 
true  that  the  required  inverse  of  each  !3j  could  be  computed  from  the  algorithm  developed 
for  the  first  order  scheme,  the  pivotal  fact  is  that  alth  ugh  Uj  is  itself  a  sparse  matrix,  its 
inverse  is  not  necessarily  sparse,  and  indeed  may  be  completely  filled  with  nonzero  elements. 
It  follows  that  the  formal  solution  of  (ol)  for  3,  given  by  the  expression 


(63)  6  =  [t  (Us-U.U^UtJ'1 

ft  (Ss-Ut  U?ct)| 

-i=i  J 

Li=i  J 

is  computationally  practical  only  when  the  dimensions  of  the  U 4  are  reasonably  small. 
Accordingly,  when  the  Ut  are  considered  to  be  indefinitely  large,  the  above  reduction 
is  untenable,  and  the  algorithm  developed  for  a  first  order  system  cannot  be  applied 
directly  to  effect  the  reduction  of  a  second  order  system.  As  we  shall  see,  however,  it 
can  be  applied  in  an  indirect  manner  to  produce  a  practical  reduction. 

Let  us  represent  Ut  more  compactly  as 


(64)  U\ 


N\ 
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FIGURE  4.  Schematic  of  coefficient  matrix  characteristic  of  a  second  order 
partitioned  system  of  normal  equations. 


FIGURE  5.  Schematic  of  coefficient  matrix  characteristic  of  a  third  order  partitioned 
system  of  normal  equations . 


As  in  Section  1,  we  let 


N»  N, 


m  -  .. 

Ln;  n, 


Mt  M, 
M\  Mj, 


By  virtue  of  the  results  of  Section  1,  a  more  detailed  representation  of  the  right  hand 
side  of  (68)  is 


Mt  M,  I"  M,  -MtQT, 

mJ  -Q,Mt  Nj*  +  Qj  MtQT, 


From  (67),  (68),  (69)  it  follows  that 


-MtQTt  UT. 

(70)  utut  u;  =  (ut  ut)  . 

t°* Mi  Ni  +Q1m1q*J  Lu\ 


which  reduces  to  V 

(71)  u,  uf  u;  =  ut  Mt  0;  -  0,  m,  q\  u;  -  u,  q-,  Mt  u\  +  ut  n"1  u\  +  ut  q,  m,  q;  u; 

This  may  be  written  as 

(72)  u1  u^u;  =  Ut  n;1  u;  +  (Ut  -  Ut  Qt)  Mt  (0,  -  Ut  Qt)T 

•“  "  —1  M 

In  a  similar  manner,  we  find  that  the  expression  Uj  Ut  c,  in  (63)  may  be  written  as 

(73)  Ut  U^c,  =  Ut  Nf  c‘j  +  (U,  -  U,  Qt)  Mt  (c,  -  Q\  ct) 

Although  it  is  not  yet  readily  apparent,  equations  (72)  and  (73)  constitute  the  key  results 
that  make  practical  the  reduction  of  a  second  order  patterned  scheme.  This  is  because, 
the  evaluation  of  certain  critical  quantities  in  these  expressions  can  be  performed  in  terms 
of  sums  of  matrices  of  low  order  that,  in  turn,  can  be  readily  evaluated  during  the  course  of 
the  reduction  of  each  of  the  first  order  systems  implicit  in  the  second  order  system.  Before 


we  demonstrate  this,  it  is  convenient  to  proceed  as  if  the  quantities  required  for  the 
evaluation  of  (72)  and  (73)  had  been  computed,  it  clarifies  matters  if  we  introduce  the 
following  set  of  auxiliaries 


(74)  [ft],  =  u,  -  u,  n;1  X)\ , 

(75)  CN 3t  =  U,  -  Ut  Qt  , 

(76)  [N]"tl=  M,  =  (N,  -  N,  N;X  NT)’1 

(77)  [a],  =s. -utN;lclf 

(78)  [c],  »6,  -QTt  c\  =c, 

As  might  be  surmised,  these  auxiliaries  have  been  formulated  to  be  analogous  to  quantities 
involved  in  a  first  order  process.  The  analog/  is  carried  further  with  the  formulation  of  the 
following  set  of  secondary  auxiliaries: 

(79)  CQ3t  =  CN3TX  CN3I 

(80)  [R]t  *=  [N],  [Q], 

(81)  [S]t  «[N]t-CR3, 

(82)  [cr]t  =  [c],  -  [Q]J  C*c]t 

As  the  [S]t  and  [c],  are  formed,  they  are  added  to  the  sum  of  their  predecessors.  This 
leads  ultimately  to  the  formation  of  the  sums 

(83)  [S]  KS]s  +  ...+CS],  , 

(84)  [c]  +[*]■  +  . ..+  • 


But  these  are  equivalent  to  the  expressions  in  bracket*  on  the  right  hand  side  of  (63), 
a  direct  consequence  of  the  fact  that 
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(85)  [S],  =Q,  -u.u,  u; , 

(86)  . 

Accordingly,  the  solution  for  6  can  be  written 

(87)  6=[M][c], 

in  which  [M}  is  defined  as 

(88)  [M]  =  [S]”1 . 


Once  6  has  been  deteimined  and  eliminated,  the  second  order  scheme  collapses  to  the 
independent  reduction  of  n  first  order  schemes.  Specifically, 


(8?)  i.-ltfs.-cco.s 


But  by  virtue  of  (62b)  and  (62c)  the  expression  Uj  c,  is  the  same  as 


the  practical  evaluation  of  which  is  precisely  the  problem  considered  in  connection  with 
first  order  partitioned  regression.  If  we  let  the  vector 


I 

I- 

t 

■  -  i 

1  .  t 

denote  the  intermediate  or  provisional  solution  obtained  by  applying  tb,e  first  order 
process  to  (90),  the  final  solution  may  be  expressed  as 

(92) 

This  shows  that  the  reduction  of  a  second  order  system  may  be  developed  in  terms  of 
independent  reductions  of  n  first  order  systems,the  solutions  of  which  are  subject  to 
subsequent  modification  once  the  vector  5  has  been  established. 

The  practicality  of  the  above  reduction  hinges  on  the  fact  thot  the  dimensions  of 
the  five  basic  auxiliaries  [N],,  [N]j,  [N],1,  [c]t  t  id  [c]t  are  determined  by  the 
orders  of  the  vectors  6  and  6t  which,  unlike  those  of  the  vectors  fij,  do  not  increase 
with  n4 .  Accordingly,  these  auxiliaries  are  not  subject  to  unlimited  growth  as  more 
and  more  data  are  processed.  However,  it  remains  to  be  shown  that  their  computation 
Is  a  practical  matter,  for  this  question  had  been  bypassed  in  the  development  beginning 
with  equation  (74). 

A  •  •  A  e 

We  begin  by  anticipating  the  result  thot  the  primary  matrices  U,f  Ut,  Nj,  ct,  e, 
that  appear  explicitly  in  the  second  order  system  (66)  are  each  expressible  as  sums  of 
matrices  generated  at  the  level  of  the  double  subscripted  matrices  Nijr  NtJ,  ctj.  Thus 
we  write 

(93a)  01  =  UU  +0„  +  ...  +0ltj  , 

(93b)  Uj  =  Ut,  +U18  +  ...  +U!ni  , 

(93c)  Nt“  +  ^18  +  •••  +  NlBi  , 

(93d)  ct  =  cu  +c13  +  ...  +c1Bj  / 

(93e)  c  t  =  c  u  +  c  l8  +  . . .  +  c ,  . 

By  virtue  of  the  partitioning  of  Uj  and  Nt  implicit  in  (62a)  and  (62b),  the  expression 

♦  *«1  T. 

N4  Uj  reduces  to  me  form 


"HOT 

\ 
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in  which  QtJ  is  defined  by 

(95)  Q.r'Nu  05, . 

If  we  then  define  the  additional  auxiliaries 

(96)  RtJ  -Djjfc,,  r 

(97)  C^],j  «014  -  ft*,  / 

it  follows  that  [N34  can  be  expressed  as 

(98)  [N],  *CN3U  +[N]l8  +  ...  +CN]Ul  r 

and  similarly  that  [c]j  can  be  expressed  as 

i 

(99)  [c],  “Cclji  +[c]1#  +  ... +[c]llt 
in  which 

(100)  [e],,  =cu  -6’,  eu  . 

Because  Ut  and  Q 4  are  row  and  column  vecto-s  of  U4  j  and  Q41,  respectively/  it  further 
follows  that  [N](  can  be  expressed  as 

(101)  [N],  =CN]U  +[N3ia  +  ...+CN]1Bl 

in  which 

(102)  CN3u=Un-UnQtJ  . 
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Finally,  we  note  that  the  matrices  [N]j  =Mt  and  [c]t  =C1  would  be  generated 
during  the  normal  course  of  the  reduction  of  a  first  order  system.  It  follows,  then, 
that  by  augmenting  the  reduction  of  a  first  order  system  by  the  generation  of 
equations  (95),  (96),  (97),  (100)  and  (102)  (all  of  which  involve  simple  operations  on 
matrices  of  low  order)  and  by  performing,  in  a  cumulative  manner,  the  summations 
(98),  (99),  (101),  one  can  generate  the  auxiliaries  (74)  through  (78)  that  are  required 
for  the  reduction  of  the  second  order  system. 

We  shall  shortly  extract  from  the  above  derivation  a  concise  computational 
outline  of  the  reduction  of  a  second  order  system.  But  before  so  doing,  we  shall 
clarify  certain  points  that  might  stili  be  obscure.  First,  let  us  return  to  the  matter 

rw  •  •  rw 

of  the  evaluation  in  (63)  of  the  expression  Ut  Ut  UJ  which  we  pointed  out  is 
impractical  to  accomplish  in  a  direct  manner  when  U  i  is  large.  The  key  to  an 
appreciation  of  precisely  why  the  Indirect  reduction  |ust  derived  is  practical,  whereas 
a  direct  reduction  is  impractical,  lies  in  a  consideration  of  the  final  term  on  the  right 
hand  side  of  equation  (72).  This  term  can  be  written  out  more  fully  as 


The  amount  of  computation  required  for  the  evaluation  of  this  expression  depends 
altogether  on  the  order  in  which  the  multiplications  are  performed.  An  inefficient 
reduction  would  consist  of  evaluating  first  the  expression  Qt  Q\  and  then  pre  and 
post-multiplying  this  result  by  Ut  and  U*.  An  efficient  reduction  would  consist  of 
evaluating  first  the  expression  D,  Qj  and  then  pre  and  post-multiplying  I^Ai  by  this 
result  and  its  transpose.  To  see  this  readily,  suppose  that  and  the  Ut }  ond  Q%i 
were  gll  scalars.  Then  one  could  easily  verify  that  evaluation  of  (1(8)  by  the 
inefficient  procedure  wouid  require  on  the  order  of  2n*  multiplications  and  n{  additions. 
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whereas  evaluation  by  the  efficient  procedure  would  require  about  2n,  multiplications 
and  nt  additions.  Accordingly,  the  practicability  of  the  reduction  of  a  second  order 
system  can  be  said  to  depend  essentially  on  efficient  manipulation  of  mafrices.  The 

~  ••  _i 

direct  evaluation  of  the  key  expression  U  t  U t  is  impractical  simply  because  it 
entails  matrix  operations  that  are  inefficient  in  view  of  the  desired  end  result  (note  here 
that  the  full  inverse  of  *Ut  is  not  generally  of  interest,  and  hence  its  evaluation  as  an 
intermediary  is  of  no  direct  value).  By  examining  the  finer  structure  of  this  expression, 
we  find  that  means  do  exist  for  its  efficient  evaluation,  and  it  is  this  fact  that  leads  to 
a  practical  reduction  of  a  second  order  system. 


It  will  be  noted  that  in  developing  the  reduction  of  the  second  order  system  we 
adopted  notation  paralleling  that  used  in  the  reduction  of  a  first  order  system  (compare 
equations  (79) - (89)  with  (50)-(57)).  Thi  suggests  that  the  reduction  of  a  second  order 
system  is  okin  to  the  double  application  of  the  reduction  of  a  First  order  system.  That 
this  is  indeed  so  is  especially  obvious  from  the  following  alternative  development 
attributable  to  John  Stephenson  (private  communication).  The  original  second  order 
system  (66)  can  be  rearranged  as  follows 
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Let  this  be  abbreviated  as  follows  in  accordance  with  the  partitioning  indicated  by 
the  broken  lines: 


P  P  5  c 

005)  -  "  .. 

P*  P  6  i 

U  J  L  J  J 

Then  the  reduced  normal  equations  in  the  vector  6  as  obtained  by  application  of  the 
first  order  algorithm  for  the  elimination  of  4"  is 

(106)  (P  -  P  P  pr)  5  =e  -T  P  c. 

It  is  readily  verified  that  the  expressions  P  P  '  7T  and  P  P  V!  reduce  to 


i  us  N?  U\  u,  nT'nj  u,  n,'‘n;  ... 

0  ...  0 

?  K** «  N8N^Ut  0  N.N^NJ...  0 


Nn  K  u; 


o  ...  n.  n^Fi; 


prlc- 


t  a,  «?*, 

s=»  _ 

*4 


n.nTc.  . 


—  ••  m1«. 

R.  N.  e. 


The  key  points  here  are  that  P  P  P  has  precisely  the  some  structure  as  P  and  that  Nt 

is  a  block  diagonal  matrix ,  The  latter  fact  permits  the  elements  of  the  above  matrices  to  be 

expressed  as  sums  of  low  order  matrices.  Thus 

O.n;1^  eQ11n^u'1, 

U,  N^Rt.  Lu,,  n;',  R;, 


(108)  N,n;1N*=  SN,,N~  N'„ 
0,  N?c,  =  SU„  c„ 

N  N^c  =  S  N  N_l  c 

,Ni  Ct  ,NI  J  ,NlJ  C!J 


—  ..-1  _ 


where  all  the  sums  run  from  j  =  1  to  n4.  The  form  of  the  reduced  normal  equations 
(106)  is 


009) 


^(Uj-d.n;1  a;) 

Ox  -Ox  Nx"1  NxT 

Ua -  U„  Ng1  UJ  ... 

O.-U.N^UT 

8 

E(8t-  UtN;lct) 

i=i  x 

Or- Ns  Ns"  u l 

Nj  -  Nj  Nj""1  NJ 

0 

fii 

1=1  _  ..  _i 

0,-NxNx  V 

UT-N8N8-l0j 

0 

N8-NaNglRj  ... 

• 

0 

6a 

• 

s 

ca-N,N^c. 

4 

• 

• 

01-n.n;1  0; 

0 

• 

• 

0 

n.-  n.n;‘ n; 

'k 

• 

• 

A  ■  a  m 

C.-N.N,  e. 

which  in  view  of  the  auxiliaries  (74)  -  (78)  may  be  written  more  concisely  os 


(110) 


tCN],  CN]x  [N]a  ...  [N]. 

A 

6 

t  ccd7 

*=i  _ 

tel 

[N]f  [Nl  0  ...  0 

[cl 

[N]J  0  [N]a  ...  0 

•  *  •  • 

«a 

• 

= 

Cc]s 

* 

♦ 

•  •  4  4 

CNJT  0  0  ...[N]n 

• 

L4-, 

• 

[c]n 

This  reduced  system  is  itself  a  first  order  system  and  its  reduction  proceeds  in  accordance 
with  equations  (79)  -  (89). 


The  above  development  makes  it  clear  that  the  reduction  of  a  second  order 
system  can  be  accomplished  be  a  double  application  of  the  reduction  of  a  first  order 
system.  All  that  remains  to  be  done  is  to  transcribe  the  reduction  into  a  concise 
computational  flow  abstracted  from  derivational  detail.  To  do  this  we  begin  with  the 
fact  that  the  system  of  observational  equations  that  ultimately  gives  rise  to  the  second 
order  system  is  of  the  form 
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>3  : 


1  I 


i  ! 


(]]})  Vu  +Blj6  +  B1)fl1  +B145n  =€tl 


i  “  1/  2/ « •  «  /  n 

i  =  •• .fftj 


The  covariance  matrix  A  associated  with  the  observational  vector  is  considered  to  be 
diagonal  ana'  partitioned  as 

(112)  A  —  diag  (Aj  As  •••  A#) 

where  in  turn 

(1 13)  At-  diag  (Ajj  Aja  •  •  •  Alnj  )• 

The  algorithm  for  the  reduction  of  a  second  order  system  then  begins  with  the  evaludrten 
of  the  following  array  of  primary  matrices 


X  X  X  cn 

fir  “ 

Bn 

(114) 

N.,  NtJ  fitJ 

= 

—ii*  •  •• 

A‘J  |Bu  Bu  Bn  <u 

y\>  NTtJ  Nti\ 

5,. 

In  terms  of  these  the  following  primary  reduced  array  is  formed 


(115) 


c*].,  CN]tJ  mti‘ 

*^il  cij 

1 

C| 

“  J 

CN]TtJ  [N]tJ  PHs, 

Nu 

R.l 

As  each  such  reduced  array  is  formed  it  is  added  to  the  sum  of  its  predecessors  thereby 
producing,  after  running  through  j  =  l,2,...,nt, 


(116) 


[N]t  [N]t  [6]t 
m\  cn]j  LcX 


S[N]1} 

£CN]tJ 


sCN]„ 


SC*]., 

s  [*•]„. 


The  provisional  solution  for  6t  is  obtained  from  this  array  by  means  of 


(117)  *>-[ N];  [C]t 

As  each  array  of  the  form  (1 1b),  is  formed  its  elements  are  operated  on  to  form  the 
secondary  reduced  array 

(118)  [cS]tC5]1)=[cN]ICi]1)-  CN]i  [N];l[[N]\C*c]J  . 

These  are  thensumme_d  as  "they  formed  to  produce  after  5  has  run  from  1  to  n,  ' 


(119)  [CS][C]|  = 


£  CS],  SCc]i 


Li=l  1=1 

A 

In  terms  of  these,  the  solution  for  6  is 

(120)  3  =  [S]"l[c]. 

The  revised  solutionsfor  each  of  the  3t  are  then  obtained  from 

(121)  =  -  [Nl^CNp,  3  i  =  l,2,...,n 

•  • 

and  in  terms  of  these  the  solutions  for  each  of  the  are  obtained  from 

(122)  «„=n;)c„-<  NJ,4t. 

With  3,  3t  and  fit,  thus  determined,  the  residuals  may  be  computed  from  the 
observational  equations  (1 1 1) 


•  •  •  • 


(123)  vt4  fti”(®uS  +  ®ijfii+®iifiii) 

Alternatively,  if  the  reduction  is  iterated  until  the  corrections  to  the  parameters  become 
insignificant,  the  residuals  may  be  computed  from 


(124)  v1J=e1J 
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where  now  denotes  the  discrepancy  vector  arising  from  the  substitution  of  the 
final  values  of  the  parameters  into  the  original  observational  equations. 

The  error  propagation  associated  with  the  adjustment  requires  the  computation 
of  those  diagonal  blocks  of  the  inverse  of  the  coefficient  matrix  that  correspond  to 
the  vectors  5,  and  *6t  j.  For  6  we  have  immediate,/ 

(125)  vcir(6)=[S]"l*[M] 

« 

For  5,,  the  appropriate  result  is 

(126)  var(6i)=:[S]71+CQ3l  [M] 

where  [Q]t  -  [N]^  [N]' .  Finally,  that  for  6, }  I,  given  by 

(127)  varC6tJ)  =  +  Q , }  var(6;)  C>\ } , 

where  Qtj  =  NtJ. 

Equations  (1 14)  through  (127)  provide  in  concise  form  the  complete  reduction  of  a 

second  order  system.  Such  a  system  can  become  indefinitely  large  in  either  or  both  of 

two  ways:  namely,  n  can  become  indefinitely  large  (as  in  a  first  order  system)  or  some  or 

all  of  the  n,  can  become  indefinitely  large.  When  the  reduction  of  a  second  order  system 

is  accomplished  by  means  of  the  above  algorithm,  the  overall  computational  effort  tends 

to  increase  linearly  with  the  dimensions  of  both  n  and  n..  Unless  the  dimensions  of  the 
•  •«  * 

Ni  and  Nij  matrices  are  grossly  disparate,  the  precise  manner  in  which  the  normal 
equations  grow  is  of  little  computational  consequence.  To  conclude  this  section  we 
would  note  that  the  three  points  listed  at  the  end  of  Section  4.1  as  characterizing  the 
reduction  of  a  first  order  system  apply  equally  well  to  a  second  order  system. 

4.4  Applications  of  Second  Order  Partitioned  Regression 

Although  the  programs  SAGA  and  GDAP  are  both  short  arc  reductions,  they  do 
not  apply  the  same  approach  to  the  solution  of  the  normal  equations.  In  GDAP,  the 
algorithm  for  a  first  order  patterned  system  is  employed.  This  is  practical  because  the 
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U,  characteristically  generated  in  GDAP  reductions  do  not  become  excessively  large. 

Thus  relatively  little  is  to  be  gained  in  GDAP  from  application  of  the  second  order 

•  • 

algorithm.  In  SAGA,  on  fie  other  hand,  the  Uj  can  become  quite  large,  particularly 
when  error  models  are  reinitialized  several  times  over  a  pass  at  many  of  the  participating 
stations.  For  this  reason,  SAGA  employs  the  algorithm  for  second  order  systems,  and  the 
program's  capabilities  are  thus  considerably  expanded  over  what  would  have  otherwise 
been  practical . 
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5.0  AUTOREGRESSIVE  FEEDBACK 


5. 1  Introduction 

SAGA  has  an  option  designed  to  take  serially  correlated  errors  into  account. 

This  capability  is  of  particular  value  when  high  sampling  rates  are  exercised,  as  in 
optical  chopping  of  passive  satellites.  Here,  successive  observations  are  likely  to 
inherit  common,  slowly  changing  errors  by  virtue  of  their  proximity  in  time  and  space. 
Even  though  such  errors  may  be  well  submerged  in  the  noise,  they  ultimately  set  limits 
to  attainable  accuracies  —  limits  that  cannot  be  overcome  by  increased  sampling. 

Thus,  we  can  assert  that  in  optical  chopping,  a  set  of  500  points  on  a  given  plate  is, 
by  virtue  of  low  order  serial  correlation,  unlikely  to  contain  significantly  more 
information  than  a  representative  sample  of  50  points.  Indeed  it  can  be  said  in  general 
that,  even  with  systematic  errors  completely  removed,  errors  in  most  channels  of  tracking 
observations  are  not  likely  to  be  strictly  Gaussian  in  character,  but,  rather,  are  likely 
to  be  serially  correlated  to  a  greater  or  lesser  degree.  In  Brown,  Bush,  Sibol  (1963) 
we  suggested  that  the  key  to  the  practical  resolution  of  such  difficulties  might  well  lie 
in  application  of  a  basic  result  of  autoregressive  theory  derived  by  Wise  in  1955  ("The 
Autocorrelation  Function  and  the  Spectral  Density  Function,"  Biometrika  42,  pp. 
151-159).  In  Brown,  Bush,  Sibol  (1964)  we  enlarged  on  this  theme  and  developed  the 
essence  of  what  we  now  call  'autoregressive  feedback.'  We  have  since  implemented  a 
limited  version  of  this  concept  in  a  successful  application  to  the  reduction  of  geodetic 
SECOR  observations.  Brown  (1966).  Not  oniy  does  autoregressive  feedback  appear  to 
provide  a  practical  answer  to  the  problem  of  serially  correlated  errors,  it  also  appears 
to  be  an  excellent  means  for  taking  into  account  the  effects  of  unmodelled  systematic 
error  which  induce  serial  correlation  into  observational  residuals  resulting  from  a 
conventional  adjustment.  We  consider  autoregressive  feedback  to  be  an  exciting 
development,  particularly  since,  as  will  be  seen,  its  implementation  can  be  so  readily- 
accomplished  as  a  natural  adjunct  to  conventional  adjustments  which  ignore  the 
presence  of  serial  correlation  Its  incorporation  into  SAGA  provides  the  program  with 
a  capability  that,  we  feel,  will  become  of  increasing  importance  as  the  significance  of 
serial  correlation  becomes  more  generally  appreciated. 


5.2  The  Autoregressive  Model 


A  stationary  sequence  of  serially  correlated  errors  ct,  Ct_lf 
is  governed  by  an  autoregressive  process  in  which 

0)  £1  “  Olj  fi«s  ***  "t  ttp  c  p  ^7 1 

where  the  a's  are  constant  coefficients  and  rj,  is  a  random  impulse  of  zero  mean  and 
variance  a3 .  According  to  this  process,  the  i  th  error  in  the  sequence  is  generated  as  a 
fixed  linear  combination  of  a  set  of  its  predr  cessors  in  combination  with  a  superimposed, 
strictly  random  impulse.  The  process  indicated  in  equation  (1)  is  said  to  be  of  pth  order 
because  p  coefficients  are  involved  in  its  description.  Specific  examples  of  the 
character  of  errors  generated  by  various  first  order  autoregressive  processes  are  to  be 
found  in  Brown,  Bush,  Sibol  £1963).  Experience  thus  far  indicates  that  autoregressive 
processes  of  low  order  provide  satisfc  :tory  stochastic  models  for  errors  encountered  in 
channels  of  tracking  observations.  In  the  case  of  Geodetic  SECOR,  sampled  at  one 
point  per  ten  seconds,  a  first  order  process  (i.e.,  ft  =at  ct„i  +  T}t)  has  been  found 
generally  to  provide  a  satisfactory  description  of  the  error  process  (Brown  (1964)).  In 
Section  5.4  we  shall  provide  statistical  criteria  for  what  constitutes  a  'satisfactory1 
autoregressive  model. 

5.3  Inverse  Covariance  Matrix  of  Autoregressive  Process 

As  we  pointed  out  in  our  earlier  work,  the  practical  utility  of  the  autoregressive 
process  in  the  adjustment  of  tracking  observations  stems  from  the  fact  that,  for  a  given 
se*  nf  autoregressive  coefficients  (which,  as  we  shall  see,  can  be  estimated  from  the 
observations),  one  can  compute  the  inverse  of  the  covariance  matrix  of  the  observational 
channel  analytically.  Equally  important,  the  inverse  turns  out  to  be  a  multidiagonal 
matrix,  the  number  of  diagonals  being  equal  to  2p+l  for  a  process  of  pth  order.  In 
Brown,  Bush,  Sibol  (1964)  we  showed  that  the  basic  result  to  this  effect  derived  originally 
by  Wise  could  be  put  into  a  more  convenient  form.  Namely,  if  A  denotes  the  covariance 
matrix  of  an  autoregressive  process  ooverned  by  coefficients  a x,  a s, . .  .,ap,  then  A  1 
can  be  expressed  as  the  following  product  of  lower  and  upper  triangular  matrices: 
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(2)  oa  A 


an-|  tfn-2  an-i  an~4 


..  oj 


-1  «1  Ofa  Ofc 

P  -1  O,  O, 

0  0  -1  «, 

0  0  0  -1 


-1  0  0  0  0 


in  which  arn-  o  for  n  >p.  If  this  result  is  applied,  for  instance,  to  o'  third  order- 
process  involving  coefficients  0f„  «*,  a, ',  one  finds,  upon  performing  the  above 
matrix  multiplication  that  A  1  assumes  the  form 


(3)  a2  A"1  - 


CU3  b_a  c_i  d  \0  0  0  o  0  0  0  0 
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o  0  0  0  0  o  \o  ...  b  db-ac-a 

0  0  0  0  0  0  0  0V-.C 

0  0  0  0  0  0  0  0  .\d 


=  1  +  a{  +  a\  o_2  -  1 1  +  tt* 

■  -«i  +  a\  «*  b-2  ■  _ai 

at  -flj 

This  illustration'  is  sufficiently  general  to  demonstrate  certain  properties  of  the  inverse 
that  are  of  pi  void  importance  to  our  concept  of  autoregressive  feedback: 

(1)  the  inverse  is  multi -diagonal,  the  number  of  diagonals  being 
equal  to  7  for  p=3  (or  equal  to  2pH  in  general) 

(2)  except  for  pxp  submatrices  comprising  the  upper  and  lower 
comers  of  the  matrix,  the  elements  of  each  diagonal  are 
constant,  the  number  of  different  constants  being  equal  to  p, 
the  order  of  the  process  , 

(3)  the  formulas  for  the  elements  of  the  matrix  are  subject  to  a 
systematic,  easily  generalized  development. 


in  which 

a  =  1  +  af  +  a}  +  a-i 

(4)  b  +  Of,  +  dfeflfe  b_, 

C  -  -Ct2  +  {*2  flfc  C_] 

d 


Expressions  for  second  and  first  order  processes  can  be  obtained  from  the  above 
development  by  successively  equating  06}  and  equal  to  zero  in  eqs.  (4). 


5.4  Estimation  of  Autoregressive  Coefficients 

Let  us  assume  that  by  some  means  one  has  obtained  a  vector  of  residuals 
(vj,  v2,  ...,  vn)  that  constitutes  an  estimate  of  the  vector  of  actual  errors 
(€),  €2,  . .  One  can  then  generate  a  set  of  autocorrelation  coefficients 
P\fP2i  •  •  •  *n  the  usua^  banner  from 

£  Vi  v*.* 

(5)  pk  - - * - 

Tv* 

*  • 
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If  we  now  write  the  autoregressive  function  in  terms  of  residuals 


(6) 


v.  = 


vt-i+Ofa  vt-a +*  ••+<*,  Vij.+Ti! 


and  regard  this  as  defining  a  set  of  observational  equations  involving  the  auto¬ 
regressive  coefficients  as  unknowns/  we  can  employ  a  least  squares  adjustment 
to  generate  a  system  of  normal  equations  for  the  determination  of  ihe  a*s. 
These  turn  out  to  assume  the  form 


(7) 


1  P,  P2  P3 

P,  V  P,  P2 

P2  P\  1  Pi 

P3  P2  P,  1 
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...  Pp_; 
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♦» 


•  • 


P  P  P  P  , 
P  p-l  p-2  p-3 


-1 

—  — 

«» 

Pi 

«2 

Pi 

Pi 

Pa 

♦ 

• 

8 

P’4 

•  , 

• 

• 

°P 

• 

PP 

!■ 

In  practice/  one  would  wish  to  establish  tne  autoregressive  process  of  lowest  order 
that  satisfactorily  models  the  observed  process.  To  do  this  one  would  begin  with  a 
first  order  process  which,  from  (7),  would  lead  to  the  estimate 

(8)  ai  35  Pi  . 

From  this  one  would  compute  the  secondary  residuals  7Jl  from 

(9)  Vt  *  €t  -  PiCi-j  . 
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If  these  secondary  residuals  turn  out  to  be  serially  uncorrelated/  the  first  order 
process  can  be  accepted  as  satisfactory.  To  determine  whether  or  not  the  Tf  s 
ore  independent,  one  would  first  compute  their  first  order  coefficient  of  corre¬ 
lation  from 

00)  .  \ 

Eul 

To  determine  whether  or  not  r  is  significantly  different  from  zero/  one  would 
employ  the  well  known  result  that  the  quantity 

(11)  In  4^ 

1-r 

Is  approximately  normally  distributed  with  mean  and  variance  of 

02)  In  ££ 

(13)  aj  =  1/(n-3) 

In  which  pis  the  true/  but  unknown/  correlation  coefficient.  It  follows  that,  at 
tho  95%  level  of  confidence/  r  will  differ  significantly  from  zero  only  if 

(14)  z  >  1.96/  Jn-3  . 

If  r  should  turn  out  to  be  significantly  different  from  zero,  one  would  then  try 
a  second  order  process.  Here,  according  to  (7),  the  coefficients  would  be  estab¬ 
lished  from 
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Pi 
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05) 


which  has  the  solution 


06) 


ai  ~  {Pi  “  Pi  Pa)/0"f?)» 
flb  «(-P?+Pa)/(WS). 


The  secondary  residuals  for  the  second  order  process  are 


07) 


<1  -  ‘"'‘.-I  *  a’!e,-2>- 


One  would  test  these  for  serial  independence  in  precisely  the  same  manner  as 
described  above  for  the  first  order  process. 

The  above  procedure  would  be  repeated  until  a  satisfactory  fit  is  obtained,  1 

■  Surprisingly,  the  simple  first  order  process  has  been  found  to  be  altogether  sufficient 
for  ranging  residuals  in  most  cases  we  have  so  far  examined.  Once  a  satisfactory 
autoregressive  model  has  been  established,  the  spectral  density  function  can  be 
computed  analytically  from  the  following  elegant  result  derived  by  Wise: 

08)  v(6)  =  i  +  0>ie'e+  tf2e‘20  +  ...  +  Q'pe‘pe){-1  +  a,e"10  +  a2e",20+ ...  + cyf1*30) 

.  where  l  -  and  8  assumes  the  discrete  sample  values 

2tt  4it  6tt  2(n-l) 


(19)  0  =  JSIL,  21. 


_  /  ™  /  ~ —  /  •  •  •/ 

n  n  n 


7T,  2tr 


-89- 


-1 
■S  1  i 

•s  l 


i 


1 


If  At  -  denotes  the  time  between  successive  points  and  T  =  n  At  denotes  the  total 
time  span,  8  may  be  put  in  the  form  »  . 

(20)  8j  *  — /  1  =  1,  2,...,n 

where  8j  may  be  said  to  correspond  to  the  frequency  of  1/2 j  At  cycle*  per 

second.  For  the  first  order  process,  eq.  (18)  reduces  to  the  following  well  known 
expression  for  the  spectral  density  function  of  the  damped  exponential  autocorre¬ 
lation  function: 

(21)  1/(6)  =  .<T2/  (1  -  2Pcos  8  +  P2) 

In  which  we  have  set  arj=P. 


5.5  Refined  Normal  Equations  (Single  Channel  Case) 

From  the  foregoing  it  can  be  seen  that  a  successful  autoregressive 
analysis  of  residuals  provides  the  solution  to  two  centrcl  problems.of  random  error 
analysis:  (1)  the  determination  of  the  inverse  covariance  matrix  of  the  observational 
vector,  (2)  the  determination  of  the  spectral  density  function  of  the  error  process. 

It  remains  to  be  shown  precisely  how  this  information  can  be  used  in  a  refinement  of 
the  adjustment  and  what  its  use  entails  in  the  way  of  additional  complication.  For 
this  purpose  we  shall  consider  the  specific  case  of  an  adjustment  of  a  single  channel 
of  observations  governed  by  a  first  order  autoregressive  process.  This  is  sufficient 
to  demonstrate  the  general  principles  of  the  operation.  Accordingly,  we  lei 

(»)  «t  «  P  «i— 1  +  V 

'  ■  '  1 

\ 
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define  a  first  order  process  in  which 


(23)  Efy)  =  0 


(24)  varfy)  =  C*  (variance  of  high  frequency  component  of  noise). 

(25)  cov  =  0  for  all  *  >  0 


Then  it  is  easily  shown  that 


(26)  E<e4) 


(27)  var  (€()  =  C*/0  “  P2)  ~  .otal  error  variance, 

(28)  «*(*,,€,_)  =  l~i, 

0*>  =  P*  • 


It  follows  that  the  covariance  matrix  A  of  on  n  vector  of  errors  is 


(30)  A  = 


pn  pn-1  pn- 2 


I'  “TMTfWIWCi*- 1  «wwi 


Let  us  now  assume  that  the  observational  vector  Is  employed  In  an  adjustment . 
Then  the  normal  equations  can  be  expressed  as  v 

(31)  (BTA"’ B)S  =  BV  c 


in  which  6  denotes  the  vector  cf  parametric  corrections  and  the  matrix  B  and  the 
vector  €  can  be  decomposed  into 


(32)  B 


where  the  b*  s  are  row  vectors  of  order  equal  to  that  of  8  and  the  e‘  s  are  scalars. 

In  a  conventional  adjustment  the  covarlanqe  matrix  is  taken  to  be  diagonal  (P=0), 
and  its  Inversion  is,  therefore,  trivial.  In  the  present  instance,  however,  the  covar¬ 
iance  matrix  is  completely  filled  (eq.(30)),  a  fact  which  introduces  complications. 
Because  A  is  in  this  instance  generated  by  a  first  order  autoregressive  process,  we 
may  employ  the  results  of  5.3  to  write  immediately  the  following  expression  for  A 


(33)  A'  * 


1 

~oT 


-  1  -P  0 

-P  .  1  +  P2  -p 

0  -p  1  +  p2 

0  0  .  -p 


0  -  ...  0 

0  ...  0 

-p  ...  o 

1+P2  ...  o 


0 


0 


0 


pippajamai  n  i  it  - 


MW* 


This  expands  to 


V 


(37)  N=-l  [(l^)(blbl+blba  +...+  bIba) 

~P  (b{  ba  +  ba  hfc  +. . .+  ba_j  ba) 

“ P  (bj  b,  +  bj  ba  +. . .+  b][  bB_i) 

-(?  (b[b,  +blba)] 

If  we  define 

(38)  Ab,  =  b,+l-  b, 
end  note  that  then 

(39)  bj_j  bt  =  b'i_j  (b1_1  +  Abj_j)  -  bt_j  b +  bj_j  Abj_j 

(40)  bl  bt_j  =  bl  (b,  -  Abt_!)  =  bl  b,  ^  bj  Abj  t  , 

It  follows  that 

(41)  b^_4  bj  +  bj  bj.j  -bj_j  bj.j  +  bj  b,  +  (bj_j  -  bj)  Abj. 

“  bi-i  bt_!  +  bj  bj  +  Abt_j  Abj_j  . 

When  this  result  is  substituted  into  (37)  and  appropriate  algebraic 
are  performed,  one  obtains 

N=  J_  [  (1-2P+P's)(bj  bt  +bj  b8  +...+  bjbB) 

aa  1 

(42)  +p  (Ab|  Abj  +  Aba  Abj +. . .+ Ab] 

+  P  d-p)(bIbj+b!bB)j 

which  is  the  same  as 

(43)  N  =  ~  [(1-P)S  BTB  +  pA  BTAB  +  p(l-p)(bjbi  +blb„)| 


x 


manipulations 


-l  Abn_j) 


-94- 


in  which 


(44)  AEH  : 


b3  “  bi 

bi,  -ba 


Ab.-J  [_bn  -  b._» j 


In  an  analogous  manner  one  obtains  for  the  right  hand  side  of  the  normal  equations 


(45)  c  =  BT  A~l  C 


=  _L  [(,-p3)  BT  €  +  PABT  AC  +p(l-p)(bl  c,  +bl  c,)] 


in  which 


(46)  '  Ac- 


Ac«_ 


*a  - 


c3  -c8 


CB  “  Ca_j 


Equations  (43)  and  (45)  are  the  results  sought.  In  analyzing  them  we  first 
note  that  when  P  =  0,  the  normal  equations  reduce  to  the  usual  expression: 

(47)  -L  (bTB)6=—  BTc  . 


Next  we  note  that  since  AB  and  Ac  are  of  first  order,  the  expressions  AB  AB  and 
AB^  Ac  are  of  second  order  and  one  can  assert  that 


I 


1 


(48) 


btb  »  abt  ab, 

BTc»ABT  Ac 


This  implies  that  as  long  os  p  is  not  too  close  to  unity,  the  normal  equations  will 

be  dominated  by  the  leading  terms  and  the  solution  vector  6  will  be  only  slightly 

dependent  on  the  value  of  p.  Hence  in  some  situations  the  solution  is  but  weakly 

affected  by  the  presence  of  even  rather  moderate  serial  correlation.  Be  this  as 

it  may,  the  covariance  matrix  of  the  solution  vector  is  very  strongly  dependent  on  the 

T 

degree  of  serial  correlation.  Even  when  B  B  dominates  the  normal  equations,  the 
covariance  matrix  of  the  parametric  vector  is  given  by 


(49)  L  & 


<1-P)S 


(BV 


Thus,  if  P  were  actually  equal  to  0.9  and  one  were  to  Ignore  this  fact  in  the 
adjustment,  the  solution  vector  itself  would  probably  not  be  very  much  affected 
(for  the  factor  (1-P)a  of  the  leading  terms  on  both  sides  of  the  normal  equations 
would  cancel  out),  but  the  covariance  matrix  associated  with  the  solution  would 
be  incorrect  by  a  factor  of  100  (or  f/(I— p)3). 

Because  of  the  relative  insensitivity  of  the  solution  vector  to  serial  correlation, 
the  residuals  from  the  adjustment  are  also  relatively  insensitive  to  serial  correlation 
(remember  we  are  considering  here  an  adjustment  restricted  to  a  single  observational 
vector).  This  means  that  it  is  a  sound  and  defensible  practice  to  estimate  the 
autoregressive  function  from  residuals  of  a  preliminary  adjustment  in  which  the 
observational  errors  are  initially  considered  to  be  uncorrelated.  The  autoregressive 
function  thus  initially  determined  can  then  be  fedback  into  the  solution  according 
♦o  the  development  of  this  section  and  the  resulting,  more  nearly  correct  normal 
equations  can  be  solved  to  generate  fresh  residuals  for  a  more  refined  autoregressive 
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analysis.  Clearly,  this  process  can  be  iterated  to  stability  (ordinarily,  a 
single  iteration  is  sufficient).  Thus  autoregressive  feedback  is  an  adaptive 
process  that  ultimately  becomes  independent  of  initial  or  a  priori  estimates  of  ’ 
noise  variance. 

Returning  to  the  normal  equations  (43),  (45),  we  note  that  the  extra 

computations  entailed  by  rigorous  consideration  of  a  first  order  autoregressive 

process  are  surprisingly  minor.  The  B^B  and  terms  would  have  to  be  computed 
T  T 

in  any  case  and  the  AB  AB  and  AB  Ae  terms,  being  of  similar  form,  follow  the 
some  logic.  Moreover,  as  is  clear. from  (42),  the  equations  can  b*s  formed  in  a 
cumulative  manner,  the  i+i  *t  step  of  which  would  involve  only  observational 
equations  from  the  i  th  and  rrj  st  observations.  Hence,  the  formation  of  the 
normal  equations  for  a  first  order  process  entails  a  scheme  in  which  a  moving  pair 
of  successive  observational  equations  are  processed  at  each  step;  similarily,  a 
second  order  process  entails  a  scheme  in  which  a  moving  triplet  of  successive 
observational  equations  are  processed  at  each  step,  and  so. on.  A  major  benefit 
to  be  derived  from  the  admission  of  an  autoregressive  process  into  the  adjustment  . 
is  the  attainment  of  more  realistic  results  from  error  propagation.  In 
particular,  it  would  permit  one  to  extract  whatever  advantages  are  to  be  gained 
from  moderately  high  sampling  rates  without  paying  the  usual  penalty  of  an  absurdly 
optimistic  estimation  of  output  accuracies. 


5.6  Normal  Equations  (Multi-Channel  Case) 


In  the  preceding  section  our  consideration  was  limited  to  an  adjustment 
.  involving  only  a  single  observational  channel.  Let  us  now  turn  to  a  multi-channel 


adjustment  in  which  each  channel  is  serially  correlated.  We  assume,  however,  that 
serial  crosscorrelation  either  does  not  exist,  or  else  is  ignored.  Thus  the  various 
observational  vectors  are  statistically  independent  of  one  another.  Let  us  now 
consider  the  problem  to  be  one  of  orbit  determination  in  which  o  set  of  error 


[ 
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coefficients  is  to  be  recovered  for  each  observational  channel.  Let  6  denote 
the  vector  ..T  orbital  parameters  and  let  denote  the  vector  of  the  error  parameters 
for  the  i  th  channel.  Let  us  assume  momentarily  that  only  the  1  th  channel  is  to  be 
exercised  in  the  adjustment.  Then  the  normal  equations  can  be  written  in  the  following 
partitioned  form 


Since  this  system  corresponds  to  but  a  single  observational  channel  we  may  apply 
the  results  of  the  preceeding  section  and  assume  that  autoregressive  feedback  has 
been  exercised  and  that  the  normal  equations  (50)  rer  le  end  product  of  the 
process.  In  this  regard  we  would  caution  that  in  practice  one  would  have  to 
exercise  temporary,  coarse  a  priori  constraints  on  the  parametric. vector  in  the  process 
of  single-channel  autoregressive  feedback  because  of  the  near  singularity 
characteristic  of  single-channel  orbital  determination  (especially  when  error  models 
coefficients  are  also  to  be  recovered).  Such  constraints  would  not  be  carried 
through  the  final  set  of  normal  equations  (50)  inasmuch  as  this  system  is  not  to  be 
solved  as  an  independent  system. 

Let  us  now  assume  that  similar,  individual  and  independent  adjustments 
exercising  autoregressive  feedback  were  performed  on  data  channels  i  °nd  * 
which  are  also  exercised  on  the  same  orbit.  The  corresponding  normal  equations 
may  thus  be  written. 


Wo  observe  that  although  equations  (50)  and  (51)  involve  different  parametric 
vectors,  the  subvector  6  is  common  to  all  three  sets.  Let  us  now  define  a  .  v 
composite  parametric  vector  for  all  three  channels  as 


We  may  then  augment  each  individual  set  of  normal  equations  with  appropriate 
sets  of  zero  elemertts  sc  that  all  will  operate  on  the  composite  parametric  vector 
Thus  we  get 
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We  arc  now  In  a  position  to  apply  the  rule  that  sets  of  normal  equations  formed 
from  independent  observational  vectors  and  involving  a  common  parametric  vector 
can  be  summed  to  produce  the  normal  equations  appropriate  to  the  simultaneous 
adjustment  of  the  merged  set  of  observations.  Accordingly  adding  (53),  (54),  (55), 
we  get 


i 


Except  for  the  fact  that  we  have  yet  to  introduce  the  a  priori  constraint  matrices 
(W,  Wt,  Wj(  Wj, )  and  supplementary  discrepancy  vectors  (e,  et,  Cj,  ek);we 
recognize  (56)  as  being  of  the  form  of  a  first  order  patterned  system.  By  following  a 
similar  development  we  could  show  that  a  second  order  patterned  system  of  normal 
equations  will  remain  unaltered  in  form  if  autoregressive  feedback  is  exercised  for 
each  of  the  observational  channels.  It  follows,  then,  that  the  exercise  of  autoregressive 
feedback  in  SAGA  does  not  seriously  Interrupt  the  general  character  of  the  data  flow. 

In  the  computational  outline  provided  in  Section  6,  we  anticipate  autoregressive  feed¬ 
back  by  the  artifice  of  proceeding  as  if  each  channel  were  governed  by  a  first  order 
process  having  a  prespecified  coefficient  of  correlation.  In  general  practice,  each 
such  coefficient  would,  for  want  of  a  better  value,  be  set  equal  to  zero  in  the  initial 
reduction.  However,  if  the  autoregressive  feedback  option  is  exercised,  the  prespecified 
values  would  automatically  be  replaced  in  a  repetition  of  the  reduction  by  values 
estimated  from  the  residuals  computed  from  the  initial  reduction  (iterated  to  convergence). 
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Although  SAGA  specifically  considers  only  a  first  order  autoregressive  process,  • 
we  would  point  out  that  in  an  article  awaiting  publication,  J.  Lynn  of  DBA  shows  that 
p  iterations  of  a  first  order  autoregressive  process  is  equivalent  to  a  pth  order  auto¬ 
regressive  process.  Thus  if  first  order  autoregressive  feedback  should  fail  to  produce 
uncorrelated  secondary  residuals,  one  could  apply  the  process  a  second  time  to  obtain 
the  effect  of  second  order  autoregressive  feedback,  and  so  on  until  serially  independent 
residuals  are  obtained. 

As  a  final  comment  we  would  point  out  that  autoregressive  feedback  is 
indifferent  to  the  origin  of  the  serial  correlation  in  the  observational  channels.  All 
or  part  of  the  correlation  accounted  for  by  the  autoregressive  process  could  consist  of 
correlation  induced  by  unmodelled  systematic  errors.  Indeed,  autoregressive  feedback 
could  in  principle  substitute  completely  for  error  modelling.  This,  however,  is 
emphatically  not  to  be  recommended,  for  deletion  of  error  coefficients  from  the  model 
would  induce  such  a  high  and  persistent  degree  of  serial  correlation  that  the  accuracy  of 
the  adjustment  would  be  severely  diluted.  Hence,  autoregressive  feedback  is  more 
properly  regarded  as  a  powerful  tool  for  accounting  for  the  combined  effects  of  natural 
serial  correlation  in  the  observations  and  induced  serial  correlation  resulting  from 
unavoidable  deficiencies  (hopefully  small)  of  the  mathematical  model. 


6.0 


computational  outline  of  saga 


6.1  Introduction 

The  analytical  basis  for  SAGA  has  been  developed  in  the  ea  Her  sections.  However, 
in  these  sections  actual  computational  steps  and  flow  are  largely  submerged  in  derivational 
detail.  The  purpose  of  this  section  is  to  extract  from  the  derivation  the  essentials  of  the 
actual  reduction.  Our  concern  here  is  with  the  heart  of  the  reduction  and  not  with 
peripheral  operations.  We  assume  all  appropriate  preprocessing  has  been  performed  for 
each  observational  channel.  Computational  details  of  various  preprocessing  routines 
available  to  SAGA  are  provided  in  Part  II  of  this  report.  Part  II  also  provides  details  of 
the  program  itself  —  set  up  procedures,  running  insl  uctions,  flow  charts,  etc.  In  addition, 
?t  outlines  results  produced  by  SAGA  for  a  tracking  network  of  twenty  stations  participating 
in  the  observation  of  twenty  seven  passes  of  GEOS  A.  The  network  exercised  the  Goddard 
laser,  the  Geodetic  Secor  system  and  various  MOTS  and  PC-1000  cameras  (Geoceiver 
observations  are  not  expected  to  become  generally  available  until  1970). 

No  attempt  has  been  made  in  the  outline  below  to  present  the  reduction  in  minute 
detail.  The  outline,  rather,  is  intended  to  serve  as  a  guide  to  a  programmer  who  in  turn 
is  guided  by  an  analyst  having  a  sound  understanding  of  the  mathematical  derivation  of 
the  reduction. 

One  point  regarding  the  outline  requires  special  clarification.  This  concerns  the 
option  for  autoregressive  feedback.  In  formulating  the  computations  we  have  found  it 
is  convenient  to  proceed  as  if  the  autoregressive  coefficients  were  known  at  the  outset.  In 
practice  this  will  not  be  the  case  and  values  of  zero  will  be  initially  employed  for  the 
coefficients.  However,  if  the  autoregressive  feedback  option  is  exercised,  nonzero  values 
of  the  coefficients  will  become  available  from  the  residuals  resulting  from  the  initial 
solution.  Thus  our  formulation  is  designed  to  anticipate  the  possibility  of  revised  auto¬ 
regressive  coefficients. 
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2a)  ory/  py  standard  schedule  number  0  for  y 

2b)  <ry,  py  standard  schedule  number  1  for  y 

2c)  orr?  py  standard  schedule  number  2  for  y 


3a)  arr  standard  schedule  number  0  for  optical  timing 

3b)  a r  standard  schedule  number  T  for  optical  timing 

2c)  <jt  standard  schedule  number  2  for  optical  timing 


4a)  <yr/  pr/  oT  standard  schedule  number  0  for  ranges 
4b)  crr,  pr,  a t  standard  schedule  number  1  for  ranges 
4c)  a,  /  P»#  <*r  standard  schedule  number  2  for  ranges 
c)  Standard  schedules  of  sigmas  of  error  coefficients. 

la)  rOctQx  standard  schedule  number  0  J 

lb)  aa/0W/ax  /  ffc/< Tt  standard  schedule  number  1  >  optical 

lc)  O&tOu'Ox*  crc,at  standard  schedule  number  2 

2a)  a,  /  0<  ,  a.  /  O',  ,  o%  ,  Ct  standard  schedule  number  0 

O  J  S  8  4  6 

2b)  0^  •  ah  •  C*a  ' °*j r  <\  ,  as  standard  schedule  number  1 

-  2c)  cr#o  ,  cr#1,  cr.3#  cr%,  cr.4,  Cte  Standard  schedule  number  2 
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d)  Standard  schedules  of  epochal  sigmas. 


1°)  OxtOrtOtt(TirOfr<Ti  standard  schedule  number  0 

lb)  a«/0rr/aSfaif0ryr9i  standard  schedule  number  1 

la)  axf(Trr<j%tOit<j;rai  standard  schedule  number  2 

«)  m  0/1/2,3,4  1  ham,onjc  coefficients  of  gravity  field 

n”  0,1, 2,3,4) 

f)  0  =  mean  rotational  rate  of  earth 
B,  Data 

a)  k  =  pass  number 

Tk  -  epoch  of  K  th  pass 

b)  X°,VJ,Zj“,Xf  ,\5/2J  =  approximate  values  of  inertial  initial  conditions  at  t”Tk 

c)  &  “  “1  /  0, 1  or  2 

=  -1  if  nonstandard  schedule  of  epochal  sigmas  to  be  used 
=  0, 1  or 2  if  standard  schedules  0, 1 ,2  are  to  be  used 

d)  If  -  "1  read  in  alternative  cr^Cy/Ci/CT^a;  #0$  for  kth  pass 

e)  :  =  ij  / i8/  i, , . . . ,i^  (mk  5  15)  schedule  of  stations  participating  on  k  th  pass 

f)  Start  and  stop  times  of  tracking  intervals  on  kth  pass  from  i  th  station  (up  to 
4  intervals  allowed). 

=  f‘rs1'  Qnd  I0**  times  of  pth  interval  (naximum  p  =4) 

g)  q  =  0  indicates  optical  data 

q  =  1  indicates  electronic  data 

if  q  =  0,  data  block  for  kth  pass  from  ith  station  consists  of  the  following. 


0  *  j  *3  jr  &  i  (time,  right  ascension,  declination)  {  =  1,2, ..  .,mt 


H)  a,  -  -1/0, 1  or  2  (note  p  =  l,2,3  or  4) 

=  -I,  If  nonstandard  schedule  of  observational  sigmas  is  to  be  used  for  pth  interval 
=  0/1  or  2  according  to  the  standard  schedule  (0, 1,2)  to  be  used  for  pth  interval 

iii)  If  a,  -  I,  use  alternative  schedule  a„,  px,  for  pth  interval 

<Jy.  P, 

M  j3p  =  -1/0,1  or  2 

=  -I,  If  nonstandard  schedule  of  error  model  sigmas  is  to  be  used  for  pth  interval 
“0,1  or  2  according  to  the  standard  schedule  (0, 1,2)  to  be  used  for  pth  interval 

v)  0p  -  -I,  use  nonstandard  schedule  era,  0^,0%,  ac ,  at 

vi)  yv  -  0  if  autoregressive  feedback  option  is  not  to  be  exercised  on  pth  interval 

=  1  otherwise 


vii)  Cj  =  nominal  focal  length  of  camera 

=  I,  data  block  for  k th  pass  from  i th  station  consists  of  the  following. 

0  tj/Tj  (time,  range)  j  *  1,2,  ...,mt 

ii)  ap,0p,yp  are  defined  as  in  optical  case  except  that  they  refer  now  to  ranging 
observations. 


C.  Computations 
Part  I.  Orbital  Integration 

Read  in  data  for  kth  pan.  Perform  orbital  integration  with  initial  conditions  given 
by  (b).  Results  of  integration  are: 


(1)  P 

(«,  »+i) 


°x  °a  •  •  •  as 

b0  ^  b,  ...  b,  l  “coefficients  of  X,Y,Z  polynomials  (degree  p) 
co  ci  c»  “•  c* 


(2)  n  = 

(a/*»(«+i)) 


in  which 


flit 

' 

ftai 


ola ...  <?l# 

ftsa  *  •  •  ft #9 

ftga  ...  fta9 


=  coefficients  of  matrizant  polynomials 


0)  -  (a,  ai  •  •  •  a,  }  =  row  vector  of  coefficients 

(1,4+1  ) 

the  matrix  Cl  can  also  be  partitioned  as 

<9,9  («+i ))  ... 

(4)  O  =  (  00  Oj  08  )♦ 

(3/#8+a)(8,as+a)(a,«s+a) 

The  matrices  P  and  Cl  are  stored  for  later  use. 


Part  I  la.  Optical  Reduction  (use  if  q  ®  0;  i  f  q  =  1  ,  go  on  to  Part  I  la) 

(5)  Use  dummy  camera  approach  to  project  each  of  up  to  four  ple  as  onto  hypothetical 
plates .  Aim  camera  axis  at  nominal  midpoint  to  get  a, a> .  Let  x  be  initially  zero,  but 
then  modify  it  so  that  final  x  axis  is  aligned  with  trace  (as  In  ;>’»nuous  trace*). 

Results  of  (5)  are  sets  of  plate  coordinates 

(6)  f  |  — 1/2,  ...,mj 


and  orientation  matrices 


A  B 
A1  B‘ 
D  E 


C 


4  =  1/2, . ..,p  (max  p  =  4) 

(p  =  number  plates  on  puss) , 


For  pth  tracking  interval  (plate)  compute  ephemeris  for  each  point  within  interval.  Thu*  for 
time  t}  compute: 


■JO 7r 


01)  n  =  o, 

j  -t, 

(»/»)'  ) 


Transform  the  matrizant  ST  to  earth  fixed  coordinates 


(12)  =  Rj  0, 

(3/e)  (3,S)  (3/9) 


From  the  Xj,Yj,Zj  obtained  from  (9)  and  the  ,Y^,Z^  obtained  from  the  master  survey 
file,  compute  for  pth  plate 


03)  \  n 


A  B  C  fXj-X, 
A'  B'  C'  .  Y,  -  V; 

D  E  F  i  l2.-^ 


and  from  these  the  plate  coordinates 


y°%,  q. 


(14)  l  no  /  =  -  (c  -  focal  length  of  camera) . 


P3  L  J  p  3 


Set  up  matrix 


fi  ^8  ^3 


fx  ft  fajp3 


'l  0  -x°Vd  A  B  C 

-1-  4  A'  B'  C' 

^p  3  0  1  -v® 0  /c 

1  n  /CJ  D  E  F 


^  f8  f. 


F  •  f*  f' 
'1  t9  rs 


AttiuKMfiNKf  a» wp  ^eagte« ; *&&& a&m t  ^ ,  -  rlivfiinr uffflii  1  i|i ii 


Using  designated  schedule  of  sigmas,  compute: 


(<Tx  +*j 

(cty  +  Y ]  Qrfc  • 


Set  up  the  Allowing  coefficient  matrix .  • 


08) 

0(D.  x 
Bp  ,(x) 

-  j:L_ 

Q  X 

(f,  f,U, 

0/0  . 

p  J 

(19) 

*■  «4  = 

cm  !  Oo 

0*0 

0/0 

l 

(j  /9) 

(lxa)  !  0xe)+i 

(20) 

«f’,w 

0/0 

1 

<*x 

p  j 

[.  1 

XP1  I 
.0/1)1 

(21) 

Bp3(x) 

1 

o . 

|  (cJ  +X3 ) 

c 

c  '  c 

Or*) 

p  1 

1 

p  3 

(22) 

fp4(x) 

.  1  i 

Ox 
*>  5 

(xoc-x  1 

{  p  J  xp  1  / 

■  • 

In  terms  of  the  above  set  up  the  composite  matrix 


(23) 


B  pjto  = 

0/1"  ^l) 


b",H  CwiOoi^M  «„b“  «  („C,M  «4>e™w 


0/»  > 


(1/3)  1(1/ «)  t(l/l)  (l/«)  (l/4) 


(*/*) 


Or*) 


•  • 


-no-  ■ 


In  which 

*  1  If  a=b 
*=  0  if  aj4>  • 

Perform  the  above  computation!  for  point  j+1  alto,  getting: 

(24)  B,/j-nW  and  c>#mW* 

Compute : 

(25)  4B,,M  ■= 

(24)  4t,,W  * 

In  terms  of  the  above  compute: 

(27)  NpJ(x)  =  (l-pMf  Bjj.WB.jW+Ps.A®;,  W^Bp1  to 

in  which  ppx  denotes  correlation  coefficient  P„  from  input  schedule.  Superscript  p  I*  added 
to  Indicate  that  in  later  solution  ps  may  differ  from  plate  to  plate.  Similarly  compute: 

* 

(28)  cP1(x)  B  0  k)«»4  to+P,»  »6<„  to. 

As  Np  3fc)  and  cp  3(x)  are  generated,  evaluate  the  sums 


b5 

(29)  N,fc)  *  £ 

!«=*» 


(30)  s(x)  -  £  r„W  +p„('*Pm)^Wc,pW4»IiW^>» 

>=*t 


6*P  (x)  and  B%,  60 

in  which  correspond  to  first  ond  last  points  on  pth  plate* 

to  and  <pt  fr) 
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N,,w  +0,. fl-P„)(B«  (x)0l,  W+Bi,  WBJ,  W) 


1 .1*  --"■■■  ■*“-  “r  ■  TfT'j-  YntfT  ftfiUtttflYrtfftt 


Computations  precisely  paralleling  those  (18)  thru  (30)  are  performed  for  the  y  coordinate. 
Thus; 


08a)  $\(y) 
(I/*) 


09a)  %\<y) 

(i,.) 


(Mb)  B f}<y) 
(*#») 


(21a)  B^j(y) 


(V  fi  « ), , 


(i,s) 


♦  t 
(*'•) 


b?;V)  i  b?;v> 

M>  I  M 


1 


C  C 


(22a)  ep  }(y)  =  -L.  {y”-/,,}. 

TPi 

In  terms  of  the  above  set  up 

Bp  ,(y)  a  same  as  (23)  with  y  replacing  x  , 

(l,XS+4p) 


Equations  (24a)  —  (30a)  are  the  same  as  (24)  —  (30)  with  y  replacing  x. 
Np(y),  c,(y).  Combine  these  with  Np(x)  and  fpfc)  to  9«t 

Np  -Np(x)  +  Np(y) 

(x8+4p)(xa-r4p) 

(31) 

c,  *  e,(x)  +  cp(y) , 

(l8+4p,  i) 
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Ultimate  result  is 


r  * 


i 


\ 


This  completes  initial  computations  for  pth  plate  from  i  th  station  for  kth  pass.  Perform 
above  for  all  plates  covering  pass  (maximum  p=4)  from  ith  station.  Add  together  all  Np 
and  Cp's.  End  result  is  Nllc ,  clk  in  which  subscripts  i  (station  number)  and  k  (pass 
number)  have  now  been  introduced.  Proceed  to  Part  III. 


Part  lib.  Electronic  Reduction  (if  cj — 1 ,  the  following  reduction  is  performed) 

For  time  t^  repeat  steps  (7)  through  (12)  of  Part  II.  Then  proceed  as  follows. 
From  computed  from  (9)  evaluate: 

(32)  r>«  *  ((X.-X'f  +(Y +(Z1-2c,f’t 

(33)  X,  -  (Xj-Jf, ti,=  VrYt)/r«“,  =  (Z.-Z'.l/r,” 

(34)  P,°  - 

(34)  e?  s‘n  ^a,eE  =  parameter  of  spheroid^ 

sin^~  Wj  =  geographic  latitude  / 

(36)  =  C(X^f  +(Yjf  +(Zcl+&Z\f$ 

(37)  ”  Xj /rj,  u |  -V^/rj,  (Zj  +  i\  Z{  )/rj  ,  direction  cosines  to  station, 

(38)  sin  Ej  =  Xj  X  j +UC!  (i  j  +i4  • 


If  ranges  are  not  corrected  for  tropospheric  refraction,  evaluate: 


(39)  f^-1 


where 


-113- 


P0  =  pressure  (mm)  at  station 

T0  —  temperature  (  K)  ot  station 

€0  ~  water  vapor  pressure  (mm)  at  station. 

For  optical  (laser)  ranges  set  e0  =0  in  (39)  and  replace  second  term  by  0.58/X3  where  Xis 
the  wavelength  of  the  light  in  microns. 


(40)  =  29.2  (T0 -30) 

(41)  a  -  2H-1)K» 

K3  =  4(H0/rcJ> 

(42)  f(E  j)  =  - - - - 

sin  E  j  +  [sin3  E  3  +  K3  ]* 

(43)  £r3  =  af(E3)  =  refraction  correction  (to  be  used  in  equation  (48)). 

If  ranges  have  been  corrected  for  tropospheric  refraction,  set  a^Oln  (43).  For  pth  tracking 
interval  set  up  the  matrices 


(44)  e“  =  ex,  u, 


oT 


V  l 


in  which  q  denotes  the  expression  (a,  +  rS  )  .  Set  up 


wi  {£■??} 


(46)  bJ}j  =  |  r  r3  r°°  rf°  f(E)jpJ  (Tt  >  fj“Tk) 

°rj> : 


(47)  eft 


(48) 


ar 


1-  {1} 


p  > 


;p  j 


-  _ ! _  /r°°  -  (r3+Arj)|  /  where  r3  is  the  measured  range. 

a'p:  1 


In  terms  of  the  above  set  up  the  composite  matrix: 


(49) 


B(l) 

i 

(1,1  7-d)) 


R(3»)  R(Sb)  r(3)  * 

BP  t  BP  1  BP  J  BP  J 


(i,a)(i,a)(i,e)(i,r>)  (i,i) 


t  RW  t  R W  t 

p  Bp  j  oas  Bp  j  -a  p  bp  j  ^4 


a(4) 

Bp  J 

(1,0 


W  t  „<*> 

p  Bp  J 


(i,i)  (i,0 


in  which 

ttb  =  1  if  a=b 
—  0  if  a/b  . 


Perform  the  .above  computations  for  point  j+1  also,  getting: 

Bp,  t+l'  *p,  ]+i  * 


Next  computations  parallel  those  Indicated  in  (25)  -  (30).  The  end  results 


are: 


(51)  Np  =  £  NpJ+ppr0-Ppp)[B^B,/P+BTb/PBWp] 

J=»p 


(^2)  Cp  X)  cp  J  +Ppr  (^"PprJC^p  f»,p  +  Bb,p  pi 


J=», 


in  which  a,p  and  b,p  denote  first  and  last  points  of  p  th  tracking  interval.  This  completes 
initial  computations  for  pth  tracking  interval  from  i  th  station.  Add  together  all  Np's  and 
cp's.  End  result  is  Nlk,  c{k  in  which  subscripts  i  (station  number)  and  k  (pass  number)  have 
now  been  introduced.  Proceed  to  Part  111. 


Part  III.  Second  Order  Partitioned  Regression 

From  Part  i la  or  Part  lib  the  matrices  Nlk  ,  clk  are  generated  from  the  observations 
of  the  k  th  pass  from  the  i  th  station.  These  matrices  can  be  partitioned  as  follows: 
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U,v  Ulk  Ulk 

(o,e)  (e,e)  (e,*,) 

Ul,  Nlk  Nlk 


(53)  Nlk  =  Ulk  Nllt  N,k 

(6,6)  (6,6)  (6^) 


u;k  n;k  Nlk+wlk 
(v)  (-tv-t) 


in  which  •{.  depends  on  the  number  of  tracking  intervals  and  the  type  of  data  (optical  or 
electronic).  Also,  the  a  priori  weight  matrix  WJk  has  been  introduced.  This  is  computed 
from: 

Plate  1  Plate  2  *  Plate  p(max  p=4) 


,  ,  ••  /  1  1  1  1  1 
(54)  WJlc  —  diag  J  2  2  2  2 

(for  gptical)  V' CX  Ot 


1111  1111 


3  3  3  3  3  3  3  3  3 

c  °b !  °x  CTc  CTw  ax  <^C 


one  for  each  unbroken 
tracking  interval 


(55)  Wlk  =  diag  , 
(for  electronic) 


fl  l±-  JL  -L  J_  J_.  JL.  J _ L_\ 

A<rJ  a*  a!  cl  a®  cr^  J 

ic y  \  *5  *a  *4  *s  »0  ‘o  *0  / 


.-(<0  ,  _  ..(0) 

5ik  -  corrections  resulting  from  q th  iterative  cycle  (6slc  =0)- 

In  terms  of  elements  of  N  k  ,  c  k  compute: 


CN]tk 

CN]U 

(e,e) 

(s,e) 

CNt 

Cnv 

je,6 ) 

(8,6) 

ulfc  uIk 


u;,  Ntk 


(Nlk+wlkr  [u;k  ■<] 


This  completes  the  separate  computations  for  the  observations  of  the  kth  pass  from  the  i  th 
station.  Repartition  the  matrices  generated  in  (56)  and  (57)  as: 


[N]tk 

j  [N]lk 

G/5) 

i  (e,6  ) 

1 

(58) 

[Kill 

:  [N]lk 

(e,e) 

!  (e,e) 

CN,,],. 

[  K  all* 

[Nils* 

(a,  a  ) 

(a,fl) 

(a,e) 

cNlai;k 

CNsa]lk 

[N8]lk 

(fl/3) 

G/  .*!  ) 

(s/a) 

[NJI* 

[Na3lk 

[N]tk 

(6/3) 

(e/3) 

(B/6) 

As  the  matrices  defined  in  (58)  are  generated  by  each  station  participating  in  the 
observation  of  the  kth  pass,  set  up  the  matrices: 


< 


(59)  [ft] 

(o«k+  3,38^4-3  ) 


[Null*  0 

0  [Nulijjk  ••• 

*  * 

m  m 

*  • 

0  0 

CNjsI^k  CNlslla* 


0 

0 

rAa3i8k 

m 

0 

[Nul,,  * 

lc 

• 

[NlA., 

"k 

CNiS]tT.k 

•k 

S  CNaali  * 
1=1  3 

_1  17- 


[NX 

(a  n.+B,  e ) 


[Nx3lak 


£  CN.V 

j=i  1 


(60)  [N1  «  [N]  +  CN1  »  +...+  CN1,  , 

(e,a)  (s,e)  (e,s)  (s,6) 


(61)  tc\ 

(a  Bk+  o,  i) 


3tj  k 

£*i3<  k 


W* 


B.fc 

1c 


i=i  J 


,  l*\ 

(«,*) 


[c3v  +C83lfk  +...+  [c]tv. 


Set  up  the  a  priori  weight  matrix  of  orbital  parameters  using  schedule  indicated  in  input: 


(62)  \  =  diag 

(s/6) 


Compute: 


-  (NX  ICNl.+wj"1  C n 


(63)  iSl  =  [NX 

(3  Bjj+3,  3  0^+3) 


r  ' 


!  8' 


m 

i  I 

f 

■  a 

!  £ 

!  ! 


(64) 


(9 


M  =  C6]x  -[N]k  {CN]X  +  Wj"1  {[c],  -Wk  £  6^!  • 

■jj+3,1)  *  -t-1  | 


This  completes  computations  for  kth  pass.  Merge  [S]k  in  the  master  survey  file 
(augmented  by  3  rows  as  columns  to  accommodate  the  coordinates  of  the  earth's  center 
of  mass.)  Repeat  the  above  process  until  all  passes  have  been  absorbed  into  master 
survey  matrix.  Final  result  is  3m+3,  3m+3  matrix  S+W  Q.V  is  augmented  to  account 
for  center  of  mass)  and  3m+3, 1  vector  c.  Compute  vector  of  corrections  to  survey  and 
center  of  mass  (superscript  (6)  is  used  to  designate  the  t  th  iteration  of  the  solution). 


(65) 


*«.) 

0 


(s  n+3  ,3m+3  ) 


1-1 

£ 

*J=% 


-  (S+W)  ^c-W  £  y  where  6^=0* 


For  each  pass,  in  turn,  compute  the  vector  of  corrections  to  the  orbital  parameters: 


(66)  = 

(V) 


[N]k  +W* 


..  ^  .(6) 
[c\  -wk  £  5k 

•  IrX 


[N3.+W, 


— _t  a  .a) 

[Nik.  [6 jk 


in  which  [6]x  denotes  that  portion  of  6  containing  the  set  of  stations  participating  on  the 
kth  pass,  i  .e. : 


6H 

* 

6. 


cU 


Joo 


The  error  parameters  for  the  i  th  station  and  kth  pass  are  computed  from: 


itiur'4  i T p'.lRlvfafcrti: ^rVAiidk' 1  *> 1 :  **ri  ^  r- 8WiVl- 4 TfnSMtK K tM  jH\VVi 


(67)  ®  =  (Nlk+W|kr|  clk  -  WIk  L  6#  -  (Nlk  +Wlk  )_1 


'IK 


’lk-  '  "  s  k  f 

M  1(00 


-t-i 

L 

t 


(0 

'  '  1  _  *  *L?- 


Add  corrections  provided  by  5^  to  survey  and  corrections  provided  by  6k  to  initial 
conditions  of  k  th  pass.  From  error  parameters  compute  the  following  corrections  for  i  th 
station  and  kth  pass. 


(68) 


6  —  crx 

X  p  *  p 

L~  i 

Mx)  E 

6^ 

P3  P3 

( 1  /  4)  4-  3- 

(V) 

Kr  l 

6  = 
yp,  yp. 

B,M  E 

3  3 

(i,4)  <6=1 

(4,i) 

6-  =  cr. 

,p3  tp3 

-l-i  1 

B,  5  lM  e 

optical  case 


In  which  the  matrices  are  based  on  the  latest  approximations  for  survey  arid  orbit.  The 
residuals  for  the  kth  tracking  interval  are  then  computed  from: 


S  *2  2 

CT.+Xj  UT 


S’} 


m  v,  » 


xj  oT 


2  ,  •  2  _a 

a  +x,  aT 


3  U*T*J 


rpj 


K-^.  +  s’} 


optical  case 


3  *2  3 


3  ff,  +  r  a 


K- vs>} 


(71) 


vtd 


2  .  ’3  2 

crr  +  r  aT 


electronic  case 
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In  (70),  (71)  the  quanlities  x^0,  y^0  and  r^0  are  based  on  the  most  recent  values  for 

J  J  $ 

survey  and  orbit. 


Compute  grand  mean  error  for  t  th  iteration  from: 


(72)  -  !  Sljm  °f  squares  of  weighted  residuals  for  oil  stations  and  all  passes 

I  degrees  of  freedom 


degrees  of  freedom  =  2(total  number  optical  points)  +  (total  number  of  ranges) 
-6(number  of  passes)  -  (total  number  error  parameters) 
-3(number  stations) . 

Solution  is  considered  to  have  converged  when: 

(73)  |aa~l)  -  aW)| 

a  prespecified  constant  or  when  five  iterations  have  been  performed,  in  the  event  the 
inequality  is  not  satisfied  for 

Part  IV.  Autoregressive  Feedback  (Option  exercised  only  if  yp  =1  for  tracking  interval) 

After  solution  has  converged,  compute  final  primary  residuals  from  (70)  or  (71). 

For  each  tracking  interval  compute: 


/  2'V*p  £vx„ 

,  =  / _ , . J -  .  , _  „  =  _ r 

v%  V  no.  pfs.  in  tracking  interval  px  no-  pfs.  in  tracking  interval 


/  Ev’»  +,vn 

I  I _ '  _  J+1  J 

'  '  \  py  yf  no.  pts.  in  tracking  interval  py  no.  pfs .  in  tracking  interval 


T>  v  r  v  . 

p  .  !  p 

> 


V  no 


.  pts.  in  tracking  interval  UPr  no.  pts.  in  tracking  interval 
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|P„  "  ur l'/iF> 

<75)  I  pB  v  =  u»//4 


Pp  r  ^Pr /  *pr 


Repeat  the  general  adjustment  with  syi(,  spy,  spr  replacing  crp>.,  afy,  cr,  r  and  with  the 
values  of  Ppx,  °py t  Ppr  from  (75)  being  used  in  place  of  values  specified  in  input.  Iterate 
to  convergence.  Recompute  primary  residuals,  as  in  (70).  In  terms  of  primary  residuals, 
compute  secondary  residuals  as  follows: 

v* p  vx.  “  Ppx  VXp 

PJ  PJ-  i 


v  =  V  -  p  V 
jv  yPj  'py  > 


vh,  = 


Compute  rms  of  primary  and  secondary  residuals  for  each  tracking  interval . 
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The  development  of  I  ho  Short  Arc  Geodetic  Adjustment  (SAGA)  program  has  provided 
satellite  geodesy  with  another  very  accureto  reduction  tool.  This  program,  wo  believe,  is  one 
of  the  moit  advanced  reduction  programs  yet  developed  for  precise  geodetic  positioning  by  mean*, 
of  satellite  observations,  It  can  handle  any  combination  of  directional  or  ranging  observations. 
Specifically,  these  are: 

a.  Optical  (aclivecnd  passive):  PC-1000,  MOTS,  BAKER-NUNN,  EC-4, 

b.  Electronic:  Geodetic  SECOti,  GEOCFIVER,  LASER, 

Each  observation  parameter  may  be  subject  to  systematic  errors  governed  by  on  error  moo'e! 
having  either  unknown  coefficients  cr  statistically  constrained  coefficients.  In  the  optical  case 
the  angular  eioivients  of  external  oricnlation  a,  a),  X  rnay  be  subjected  to  slight  adjustments  that 
are  consistent  with  their  actual  accuracies.  The  orientation  accuracy  is  typically  0 '.'5  of  arc  in 
a,  u>  and  2'1Q  in  k ,  The  exercise  of  orientation  error  model  constraints  is  particularly  important 
in  chopping  shutter  observations.  A  conventional  reduction  of  a  large  quantity  of  chopping  shutter 
observations  par  plate,  with  the  existence  of  systematic  error  in  orientation  would  lead  to  cn  unduly 
optimistic  error  propagation.  In  a  chopping  operation  (as  opposed  to  a  flashing  light),  allowance; 
should  be  made  for  uncertainties  in  in  tor-station  liming  synchronizut'on.  The  inter-station  timing 
bias  simulates  the  case  of  non-synchronization  between  stations  observing  a  common  satellite  pass. 
The  electronic  error  model  accomodates  biases  in  zero  sot,  timing,  frequency  (satellite  oscilblcr), 
frequency  (ground  oscillotor),  frequency  drift  and  refraction,  SAGA  will  accept  as  many  as  4 
observation  intervals  from  the  various  trackers.  This  allows  observation  drop  out  with  a  new  zero 
setting  of  rouges  and  re-orientation  of  the  optical  tracker. 


*  The  material  presented  in  the  first  two  sections  of  Part  II  was  originally  given  in  a  paper  at 
the  American  Geophysical  Union  (AGU)  Meeting  in  Washington,  D.C.,  in  May  1969  as  a 
joint  effort  by  DBA  and  AFCRL  (see  Reference  6). 
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The  error  model  cdopted  foi  ranging  syslemt  is  of  the  form 

6r  -  OjiQjH  0^ T:1  +  a,,r  -t  a4f  -<  aG  esc  E 

in  wliich 

r,  r  -  range  and  range  rate  at  time  rlr-Oat  epoch) 

E  =  iocal  elevation  angle 

ao'aJ.r aS' aa/an  ae  -  error  coefficients  accounting  for  systematic  errors  in  zero 

set,  frequency  offset  (between  satellite  and  ground  station 
oscillators),  frequency  drifts,  frequency  bias,  timing  bias 
and  residual  refraction. 

The  model  was  derived  to  apply  specifically  to  Geocciver  observations  but  is  applicable  to 
ranging  systems  in  general  when  appropriate  constraints  are  placed  on  the  error  coefficients. 

Systematic  eriors  in  opticai  observations  ore  assumed  to  be  governed  by  error 
models  of  the  form: 


&x  ~  Oo  fi  *  osfa  +  a3fa  r  a4f4  +  aGf6 

6y  =  <bfi  +  aafi  +  aaf^  +  a4f4  +  apf^ 

where 

-  “(c3  +x3)/c  =  xy/ c 

fs  “  xy/c  ,s  -  c 

=  y  f'3  =  -x 

h  =  x/c  f;  =  y/c 

=  «*  =  y 

in  which 

x,y  ~  plate  coordinates  of  satellite  image 
x,y  =  rate  of  chonge  of  plate  coordinates 
c  —  nominal  focal  length  of  camera. 
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.’...  v’  vOvVoii'.viloi  01  c  generally  reconstructed  from  published  angle:,  (right  ascension, 
uecdn^iioi.)  by  a  process  of  ’dummy  camera  projection1,  wherein  the  given  angles  arc 
projected  or, to  the  plate  of  a  fictitious  camera  aimed  at  the  middle  of  the  recorded  arc. 
The  optical  e i for  coefficients  a0  through  a.;  account  for  the  combined  effects  of  bioses 
in  tile  angular  elements  of  orientation  of  the  camera,  the  elements  of  interior  orientation 
(coordinates  of  principal  point  and  focal  length),  and  timing.  They  also  account  for  any 
linear  drift  in  the  direction  of  the  camera  axis  throughout  the  exposure. 

To  demonstrate  the  capabilities  and  the  performance  of  SAGA,  we  reduced  a 
network  consisting  of  optical  and  ranging  observations.  These  involved  the  PC-1 000 
(active  and  passive),  MOTS,  SECOR  and  the  LASER.  Twenty  seven  (27)  orbirs  were 
chosen  that  provided  a  good  geometrical  relationship  between  observing  stations  and 
observed  orbits  for  accurate  station  position  determination  (see  Figure  1). 
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Twenty  stations  pvn  tici|xjted  in  the  complete  network.  The  one  sigma  error  of  the 
network  survey  determination  was  typically  3  to  5  meters,  depending  upon  the  station 
geometry  relative  to  the  network  and  the  amount  of  observations  from  the  station.  The 
error  propagation  is  internal  to  the  network  with  the  origin  at  Hunter  AFB. 

The  optical  obsen  ationr.  were  assumed  subject  to  orientation  biases  of  1"  0  in  a,  U> 
and  2'.'0  in  X  for  all  posses.  The  orbit  parameters  were  completely  relaxed  to  1x10*  meters 
in  position,  5,  mcters/sccond  in  velocity.  The  range  measurements  were  assumed  to  be 
subject  to  a  systematic  bias  of  10  meters  with  a  random  noise  of  1  meter.  The  random  error 
of  the  optical  observation  was  6  microns  in  plate  x  and  y  coordinates.  See  Tablo  I  for  all 
input  errors. 

The  SECOR  observations  were  made  on  orbits  that  had  no  optical  observations.  Four 
orbits  were  observed  by  four  SCCOR  stations  (5333  ,  5001,  5649,  5861)  and  one  orbit  was 
observed  by  three  stations  (5333,  5001,  5649).  Matching  of  available  SECOR  observations 
with  available  optical  measurement .  resulted  in  little  strength  from  either.  Therefore  the 
five  were  selected  that  provided  good  geometry  with  the  four  SECOR  stations.  The  SECOR 
stations  were  incorporated  into  the  overall  network  by  applying  constraints  between  the 
SECOR  and  optical  co-Iocaicd  stations  (e.g.  3861,  5861,  see  Figure  2). 

Laser  range  measurements  from  station  7051  were  obtained  on  four  passes  in  conjunction 
with  optical  observations.  Optical  station  1042  and  station  7051  are  again  treated  as  co¬ 
located  stations  by  constraining  the  direction anddistance  between  them. 

Table  II  reflects  the  strength  of  each  pass  into  the  overall  solution.  For  instance, 
orbits  twenty  one  and  twenty  two  are  very  weak  in  orbit  determination  since  only  two  stations 
observed.  The  accuracy  of  the  orbit  determination  reflects  the  degree  of  contribution  of  the 
pass  to  overall  surveys.  The  accuracy  of  the  orbit  position  (X,  Y, 7.)  is  largely  dependent  on 
the  intersection  angle  of  station  observations  and  the  accuracy  of  the  measurement*.  The 
determination  of  velocity  components  improves  significantly  as  the  length  of  the  data  span 
observed  increases. 


TABLE  1 


Station  Participation 


Orbit  1N0. 


Participating  Stations 


Intervals 


a 


Observed  Sets 


9 

1034, 

7039, 

7037 

2 

1042, 

7037, 

7040 

35 

7039, 

7040, 

1042 

5 

7040, 

7036, 

7039 

36 

7036, 

1022, 

7075 

37 

7039, 

7037, 

7075 

1 

7040, 

7037, 

7039 

3 

7040, 

7037, 

1034 

8 

7036, 

7039, 

7037 

10 

7040, 

1042, 

7039 

11 

7039, 

3405, 

3402 

12 

7039, 

1022, 

7040 

13 

7036, 

7037, 

7039 

4 

7037, 

7036, 

7039 

20 

3657, 

3861, 

3401 

21 

3401, 

3405 

22 

3401  > 

3106 

23 

3657, 

3405, 

3648 

24 

3657, 

3106, 

3648 

25 

3657, 

3405, 

3402 

26 

3401, 

3402, 

3648 

27 

3402, 

3401, 

3657 

40 

5001, 

5649, 

5333 

29 

5001, 

5649, 

5333 

42 

5001, 

5649, 

5333 

41 

5001, 

5649, 

5333 

43 

5001, 

5649, 

5333 

7075  4  9 

7039,  1022,  7075  3  12 

7051  2  6 

1022,  7051  2  6 

7051  2  4 

7040,  7051  4  9 

2  4 

7036  ,  7039  ,  7075  3  11 

1022,  1042,  7075  3  9 

3  6 

3861  2  5 

1042,  3861  3  8 

2  4 

7075  4  10 

1  3 

2  2 

2  2 

2  5 

2  4 

1  3 

3106  2  5 

3106  2  4 

5861  1  4 

5861  1  4 

5861  1  4 

1  3 

5861  1  4 


Total  number  of  observation  sets  131 


*The  number  of  observed  sets  represent  the  total  number  of  camera  photos  and  intervals 
of  range  observation  for  each  pass.  Each  flash  sequence  represents  one  interval. 


TABLE  2 


A  Priori  Constraints  on  Error  Parameters 


Station  No. 

Latitude 

Input  Sigmas 
Longitudo 

Height 

(sec.) 

(sec.) 

(maters) 

5001 

0.200 

0.20G 

5.00 

5330 

0.200 

0.200 

5.00 

5,-y'O 

0.001 

0.001 

0.10 

0OO  1 

0.200 

0.200 

5.00 

7051 

0.200 

0.200 

5.00 

3405 

0.300 

0.300 

8.00 

3  1-02 

0.200 

0.200 

5.00 

3557 

0.200 

0.200 

5.00 

3106 

0.350 

0.350 

8.00 

5061 

0.200 

0.200 

5.00 

3401 

0.200 

0.200 

5.00 

7040 

0.350 

0.350 

8.00 

1022 

0.200 

0.200 

5.00 

1054 

0.200 

0.200 

5.00 

1C  42 

0.200 

0.200 

5.00 

7037 

0.200 

0.200 

5.00 

7  036 

0.200 

0.200 

5.00 

7039 

0.800 

0.800 

25.00 

7075 

0.200 

0.200 

5.00 

3643 

0.001 

0.001 

0.10 

Observations 


Optical: 

Orientation: 

aa  -  1 .  Arc  seconds 

U 

=  2' 

Timing  (inter-station): 

_  .  1/r.4  (aetive) 

Ct  =  1  x  lU* 

Ot  =  3  milliseconds  (passive) 

Measurement: 


ox  -  5  microns  \ 

,  .  >  Plate  coordinat 

ar  -  5  microns  J 

q.  10  microns  Poca!  length 

at  *=  1  x  10'4  Random  timing 


Or  -  5  mien 


Plate  coordinates 


Electronic; 

Zero  set,  a,  =  10  meters 
Random  Range,  ar=  1  meter 
Timing,  at~  1  x  10-4 


Position,  Velocity 

~  Qy  =  O7  1  x  10  4  meter 
at  ~  oj-  oj  meters/seconds 


TABLE  3 

Standard  Deviation  of  Recovered  Orbital  Elements 
(Meters  and  Meters/Second) 


Pass  No. 

Px 

az 

0^ 

9 

14 

14 

37 

.027 

.031“ 

TosT 

2 

9 

12 

7 

.029 

.033 

.026 

35 

8 

9 

8 

.02) 

.024 

.02) 

5 

34 

28 

34 

.032 

.026 

.040 

36 

6 

8 

7 

.031 

.022 

.020 

37 

7 

8 

7 

.024 

.023 

.027 

« 

i 

69 

272 

279 

.158 

.411 

.468 

3 

11 

14 

8 

.022 

.028 

.021 

8 

12 

19 

14 

.035 

.054 

.040 

10 

15 

17 

19 

.035 

.038 

.064 

11 

27 

30 

12 

.067 

.123 

.132 

12 

14 

7 

8 

.037 

.036 

.044 

13 

29 

24 

44 

.088 

.103 

.101 

4 

10 

17 

10 

.019 

.031 

.021 

20 

7 

13 

13 

.518 

.887 

.944 

21 

29 

92 

184 

.397 

.394 

1.097 

22 

128 

140 

223 

.786 

.936 

.906 

23 

20 

9 

19 

.034 

.038 

.043 

24 

78 

23 

20 

.347 

.044 

.140 

25 

15 

18 

23 

.833 

1.126 

1.634 

26 

11 

15 

13 

.050 

.083 

.115 

27 

5 

11 

8 

.196 

.321 

.432 

40 

29 

31 

12 

.05? 

.094 

.052 

39 

30 

19 

18 

.111 

.078 

.062 

42 

46 

10 

22 

.048 

.036 

.029 

41 

25 

54 

10 

.061 

.105 

.070 

43 

17 

36 

26 

.037  . 

.035 

.054 

*"8- 


TABLE  4 


Station  Coordinate  Adjustments  (Meters) 


Stut ion 

Corrections 

Standard  Error 

AX 

AY 

AZ 

5001 

2.59 

-3.38 

-5.03 

3,4 

4.5 

3.2 

5333 

5.77 

-0.40 

-12. \b 

3.9 

4.8 

4.7 

5649 

0.0 

0.0 

0.0 

0.03 

0.09 

0.06 

5S61 

8.51 

-8.50 

-5.02 

3.3 

3.3 

3.2 

7051 

-7.99 

-9.39 

-8.11 

4.1 

3.v 

3.5 

3405 

-2.0 

-9.66 

-15.32 

7.2 

7.4 

7.9 

3402 

10.41 

-0.42 

-2.77 

4.9 

4.9 

5.4 

3657 

-3.02 

6.29 

-4.93 

3.1 

4.3 

3.2 

3106 

15.57 

9.76 

5.00 

9.2 

8.0 

9.7 

3361 

8.48 

-8.56 

-5.00 

3.3 

3,3 

3.2 

3401 

-7.34 

-3.98 

4.37 

4.0 

5.0 

4.8 

7040 

-2.42 

23.81 

-13.91 

6.2 

6.5 

6.9 

1022 

-1.16 

2.03 

-2.72 

•  4.4 

4.3 

4.8 

1034 

-1.81 

-6.31 

0.69 

4.0 

5.0 

5.2 

1042 

2.94 

2.34 

2.74 

4.0 

4.0 

3.5 

7037 

-2.28 

-0.68 

-10.06 

3.7 

4.1 

4.1 

7036 

13.77 

-13.58 

3.58 

4.7 

4.6 

4.7 

7039 

18.71 

5.74 

6.81 

5.5 

6.0 

5.5 

7075 

-0.60 

5.53 

2.07 

3.5 

4.4 

4.3 

3648 

0.0 

0.0 

0.0 

0.03 

0.09 

0.06 

3657 


mm 


The  totai  corrections  (Table  IV)  to  survey  and  the  propagated  errors  are  consistent 
with  the  observation  network  with  the  origin  at  Hunter  AFB.  The  propagated  errors  reflect 
the  internal  network  accuracies  and  illustrate  the  strength  of  this  particular  solution.  The 
recovery  of  some  stations  was  weaker  than  others  due  to  the  limited  observations  from  the 
station  Bermuda  was  a  priority  factor  in  selecting  orbits,  therefore  the  network  provided 
strong  geometry  for  the  determine  on  of  its  coordinates.  Bermuda  (7039)  participated  in 
thirteen  (13)  passes  which  «nvo'/  <d  twenty  four  (24)  plates  (or  24  different  flash  sequences). 
Antigua  participated  in  only  four  different  passes  with  a  total  of  5  plates.  The  position 
recovery  tor  Antigua  was  nearly  twice  that  for  Bermuda.  On  the  other  hand  the  input  error 
for  Antigua  was  much  smaller  than  that  for  Bermuda  (meaning  the  a  priori  information  was 
better).  The  improvement  of  the  recovered  coordinates  over  the  a  priori  information  was 
not  as  significant  for  Antigua.  Table  III  lists  the  observing  station  for  each  pass. 

A  typical  recovered  bias  in  optica!  orientation  a,  to,  and  x  was  .2  arc  seconds. 

The  largest  recovered  bias  was  .8  arc  seconds.  Recovered  biases  in  laser  observations  ranged 
from  4  to  10  meters  and  SECOR  biases  ranged  from  6  to  27  meters.  The  random  noise  was 
approximately  1  to  2  meters  in  both  LASER  and  SECOR  measurements. 

This  program  accommodates  adjustments  to  the  center  of  mass.  The  results  presented 
were  obtained  with  the  center  of  mass  held  fixed.  As  a  preliminary  experiment  an  adjustment 
was  mode  with  the  center  of  mass  relaxed  to  50  meters  in  X,  Y,Z.  The  corrections  to  the 
coordinates  of  the  center  of  mass  were  10,  35  and  -27  in  X,Y  and  Z  respectively  with 
standard  devicitiens  of  39  meters.  The  change  in  station  surveys  was  negligible.  This 
experiment  was  made  to  demonstrate  the  programs  capabilities  in  recovering  the  center  of 
mass.  The  network  is  concentrated  on  a  small  section  of  the  earths  surface  and  the  geometry 
with  respect  to  the  center  of  mass  was  very  weak.  Further  studies  will  be  undertaken  with 
emphasis  placed  on  orbital  geometry  and  observing  station  locations  relative  to  the  center  cf 
mass. 
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conclusions 


1:1  view  of  t lie  various  instrumentations  and  observing  modes  used  in  this  reduction, 
it  's  apparent  that  tins  program  could  bo  used  with  any  type  optical  or  ranging  system 
presently  employed  in  global  tracking.  The  Gcoceivcr  was  given  some  special  consideration 
in  the  error  model  of  range  measuring  type  systems.  Terms  for  frequency  biases  and  frequency 
drift  were  designed  to  accommodate  errors  introduced  by  the  satellite  and  ground  oscillators. 

This  presentation  is  intended  to  primarily  demonstrate  the  capabilities  of  the  SAGA 
reduction  program.  It  is  not  our  intention  to  imply  that  the  results  should  lead  to  a  change 
in  position  of  Bermuda  or  any  other  station  in  the  network.  But  the  corrections  and  accuracies 
in  the  subject  solution  arc  meaningful  with  respect  to  this  network.  The  accuracy  of  5  to  6 
meters  recovery  of  the  Bermuda  coordinates  clearly  illustrates  the  potential  of  recovering 
the  surveys  of  other  stations  that  have  large  survey  errors.  The  reduction  demonstrates  the 
impressive  potential  of  SAGA  as  a  tool  for  establishing  continental  and  inter-continental 
surveys  to  a  high  degree  of  accuracy. 
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The  SAGA  program  is  a  multi -orbital  reduction  program  designed  to  recovery  station 
geodetic  cooibinaies  and  bias  parameters  associated  with  the  measuring  system  being  employed. 
As  many  as  four  zero  settings  of  the  defined  orbit  may  be  used.  This  provides  the  capability 
of  observing  as  many  a;  four  different  segment  combinations  of  the  orbit  by  the  various  trackers. 
Optical  and  ranging  observations  are  accepted  and  these  are  specifically  PC-1000,  BC-4, 
MOTS,  BAKER-NUNNS,  SECOR,  LASER,  and  GEOCEIVER. 

Input  formats  are  primarily  the  same  as  the  GEOS  format  except  in  special  cases  in 
which  the  data  has  been  processed  through  a  data  prep  program  (see  Appendix  D)  for  special 
corrections.  This  is  a  program  option  which  was  included  to  provide  program  simplicity  and 
more  importantly  to  provide  SAGA  with  consistent  and  pre-edited  data. 

SAGA  consists  of  six  major  programs  which  are  directed  by  a  control  program.  The 
first  two  (MASTER  and  PREP)  read  all  input  data  and  store  them  into  disk  files  or  tape  for  use 
in  the  iteration  cycle.  They  also  compute  the  station  covariance  matrix  and  baseline 
constraints.  The  orbit  integrator,  which  is  part  of  the  iteration  cycle,  performs  orbit 
integration  from  the  time  of  the  initial  conditions  to  the  desired  epoch.  At  the  time  of 
epoch  it  updates  the  position,  velocity  and  time  on  the  first  iteration.  Subsequent  iterations 
need  only  to  be  expanded  on  the  corrected  orbit  and  no  further  integration  performed.  The 
user  should  provide  initial  conditions  as  near  the  desired  epocn  as  possible  to  prevent  error 
build  up  which  leads  to  additional  iterations  to  acquire  convergence.  NORMAL  forms  the 
normal  equations.  SOLVE  solves  for  survey  correction.  OUTER  solves  for  orbit  and 
observation  biases  arid  updates  the  master  file. 

Five  iterations  arc  the  maximum  number  of  iterations  before  cutoff.  If  convergence  is 
not  obtained  before  the  cutoff  point  the  data  is  probably  unstable  and  may  be  diverging.  At 
thi.  time  one  should  analyze  the  data  characteristics. 
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4.0  SUBROUTINE  INSCRIPTION 


This  section  gives  a  brief  description  of  subroutines  and  their  primary  mathematical 
operation.  The  calling  sequence  of  each  is  defined  as  follow' 


This  subroutine  clears  an  array  X  of  size  N. 

X  Array  to  bo  cleared 

N  Number  of  elements  in  X 

AUGVRT  (A,  NR,  NC.  NO) 

AUGVRT  replaces  the  matrix  A  by  A  inverse  augmented  by  the  solution  matrix. 

A  The  matrix  to  be  inverted 
NR  Row  dimension  of  A  to  be  inverted 
NC  Column  dimension  of  A  to  be  inverted 
ND  Actual  dimension  of  A 

MATMPY  (A,  NRA,  NCA,  B,  NRB,  NCB,  C,M1  ,M2) 

MATMPY  will  multiply  A  and  3  and  store  into  C.  The  options  of  Ml  provide  any 
allowable  transpose  combination  and  Ml  allows  such  storage  combinations  as  the  simple 
product,  negative  product,  sum  into  previous  C  (positive  or  negative). 

A  First  array  to  be  multiplied 
NRA  Rows  of  A 
RCA  Columns  of  A 
B  Second  array 

NRB  Rows  of  B 
NCB  Columns  of  B 
C  Storage  array  of  product  AB 
Ml  Transpose  option 

M2  Action  of  product  R  fo  storage  array  C 
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UPDATE  (iCN,  KTR,  I  PS,  IX I  ,  XPO,  YPQ,  ZPO,XOT,  YOT,  ZOT) 
This  program  updates  the  matrizant  for  a  new  step  in  time. 


ICN 

Makes  decision  on  to  integraic  posilion  and  velocity 

KTR 

Number  of  turns  to  use  in  series 

EPS 

Truncation  error  limit 

DEL 

Valid  step  size 

XPO,YPO,ZPO 

Ir.put  position 

XOT,  YOT,ZOT 

Output  position 

ZER  (A,  N,M) 

ZER  zeros  the  array  A  for  N  times  M  values  starting  at  location  A(l). 

A  Array  to  be  set  to  zero 

M  Rows  of  A 

N  Columns  of  A 

ATRACT  (L,  K, B, S,  IQT,  NS,  NROW) 

This  routine  will  disect  array  B  in  3x3  matrices  C  (3x3)  and  write  a  file  on  tape  or 
disk  identified  by  station  ID  from  array  L. 

L  Array  containing  identification  of  disected  3x3  arrays 

K  Number  of  3x3  matrices  along  diagonal  of  B 

B  Array  containing  K  squared  3x3  matrices 

S  Array  containing  K  3x1  matrices 

10T  Unit  number  of  tape  or  disk  to  be  written  on 
NS  Dummy  integer 

NROW  Rank  of  desired  partitioned  size  array 

DUMMY  ((?  aST,OR,A,E,TP,IC,IS,PH,  CH,  HE,IDNT,P,Q) 

This  program  simulates  a  dummy  camera  projection  of  right  ascension  and 
declinations  to  plate  coordinates. 
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GAST  Greenwich  sidereal  time 


OR  Output  orientation  matrix 

A  Array  of  plate  X  coordinates 

E  Array  of  plate  Y  coordinates 

TP  Time  of  i  th  X,  Y  coordinate 

IC  Number  of  coordinates  read 

IS  Station  identification  number 

PH  Array  of  station  latitudes 

CH  Array  of  station  longitude 

HE  Array  of  station  heights 

IDNT  Array  of  station  identification  numbers 

P  Latitude  of  observing  station 

Q  Longitude  of  observing  station 

COVAR  {".A,  SLAM,  SH,  R,  A,  ECC,  SIGMA) 

This  routine  computes  the  covariance  matrix  of  a  station  in  geocentric  coordinates 
given  geodetic  error  in  latitude,  longitude  and  height. 

SLA  Latitude  of  station 
SLAM  Longitude  of  station 
SH  Height  of  station 
R  Output  covariant  matrix  (3,3) 

A  Semi-major  axis  of  earth 

ECC  Eccentricity  of  earth 

SIGMA  Input  errors  in  latitude,  longitude  and  height 
EXTRAC  (L,  K,  B,CC,  NS) 

This  subroutine  extracts  3x3  matrices  from  array  B  and  wrifes  them  on  unit  3. 

L  Array  containing  position  of  3x3  matrices  in  B  that  are  going  on  unit  3 

K  Number  of  3x3  matrices  along  the  diagonal  of  B 


3  Array  being  disectcd 

CC  State  vector  to  be  disectcd  into  3x1  matrices 
NS  Size  of  submatrix . 

DRIVER  (NN/1U,TM,GH,DT/XQT,VEM) 

Th.s  program  evaluates  the  coefficients  for  position,  velocity  and  matrizant  at 
time  TM. 

NN  If  NN  =  1  read  tape  unit  IU  for  coefficients  and  orbit  parameters  at 
epoch,  otherwise  expand  only  about  epoch 

1U  Tape  unit  be  read 

TM  Time  of  observation 

GH  Corrected  output  rime  from  epoch 

DT  Increment  in  time  that  coefficients  are  good 

XOT  Output  position  and  velocity 

VEM  Output  matrizant 

MATRUP  (KTR,  DEL,UVM,UVQ) 

Subroutine  to  update  matrizant  with  respect  ro  time. 

KTR  Number  of  terms  in  series 
DEL  Increment  in  time  from  epoch 

UVM  Coefficient  matrix 

UVO  Updated  matrizant 

CRV  (AX,  ECC,  VC,  IDNT,  EPS,  ID,  I  PR) 

CRV  updates  station  geocentric  coordinates  from  corrections  in  previous  iteration.. 


h, ' 

Semi-major  axis  of  earth 

ECC 

Eccentricity  oLecutb 

VC 

'Array  containing  station  X,Y, Z 

IDNT 

Array  of  identification  numbers 
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EPS  Array  of  station  corrections 

iD  !D  of  station  match 

I  PR  Fiog  set  if  station  has  already  been  outputed 

LLH  (AX/IISQR,P/Q/R/OP/OL/OH) 

jeocentric  Cartesian  coordinates  to  geodetic  coordinates. 

Scmi-rnajor  axis  of  earth 
Eccentricity  squared  of  earth 
Geocentric  Z 
Geocentric  X 
Geocentric  Y 
Output  latitude 
Output  longilude 
Output  height 

ATANN  (X,Y) 

This  is  a  function  that  computes  the  angle  between  the  vectors  X  and  Y. 

X  First  vector 

Y  Second  vector 

DEG  (AN,1,J,S) 

DEG  converts  an  angle  AN  to  degrees/ninutes  and  seconds. 

AN  The  angle  in  radians 

1  Output  degrees 

J  Output  minutes 

S  Output  seconds 


Conv 

AX 

ESQR 

P 

Q 

R 

OP 

OL 

OH 
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GEOC  (XL/YL<  HT,  A,ES,V) 


GEOC  convcrfs  geodetic  position  ('I',X,h)  to  geocentric  coordinates  (X,Y,Z). 

XL  Latitude  (radians) 

YL  Longitude  (radians) 

HT  Height  (meters) 

A  Semi-major  axis  of  earth  (meters) 

ES  Eccentricity  of  earth  squared 

V  Returned  array  containing  X,Y,Z 

ERCOE  (EPP,  NIP,  K,  L,  J) 

ERCOE  outputs  observation  bias  corrections  of  electronic  error  model  terms. 

EPP  Array  containing  bias  corrections 

NID  Array  containing  station  identifications 

K  Defines  observation  interval 

L  Station  identification 

J  Index  defining  station  in  pass 

ERMOD  (EPP,  NID,  K, L,  J) 

ERMOD  outputs  observation  bias  corrections  of  optical  error  model  terms. 

EPP  Array  containing  bias  corrections 

NiD  Array  containing  station  identifications 

K  Defines  observation  interval 

L  Station  identification 

J  Index  defining  station  in  pass 


V 


XPO 
YPO 
ZPO  ] 

cnmI 

SNM  j 

LOT 

iCT  I 
i 

UMT 

VMT 

CTB  ( 
CTT  j 

ERD 

XMU 

ALF 

OMG 

ECC 

NTE 

KTR 

KDR 

NHT 

CDC 

CTW 

KEY 

DMT 

KRG 


Power  series  for  position  (X,Y,Z) 

Potential  coefficient 

Control  tables 

Expansion  series 
Expansion  series 

Coefficient  tables 

Radius  of  earth 
Gravity  constant 
Greenwich  hour  angle 
Rotation  rate  of  earth 
Eccentricity  squared 
Tables  control  parameter 
Number  of  terms  in  series 
Number  of  terms  in  drag 
Zero 

Ballistic  coefficient 
Drag  time 

Integration  control  constant 
Dummy  array 

Defines  the  number  of  terms  in  series. 
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A.  Input 


C 
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This  section  will  describe  the  program  input  parameters  and  illustrate  the 
output  formats.  Input  will  be  defined  as  specification  cards  and  a  complete 
orbit  sot.  A  run  consisting  of  more  than  one  orbit  is  made  by  simply  stating 
the  number  of  orbits  on  the  first  specification  card.  A  complete  set  of  N  orbits 
is  illustrated  in  Figure  1 . 


Baseline  Constraints 


Card  No. 

Name 

Field 

Columns 

Units 

Description 

1 

IP 

14 

1-  4 

Station  ID  (1st  station) 

IQ 

14 

5-  8 

Station  ID  (2nd  station) 

IS 

12 

9-10 

7^0,  read  linear  constraints 

AL 

F9.0 

11-19 

Degrees 

A  priori  azimuth 

EL 

F9.0 

20-28 

Degrees 

A  priori  elevation 

RA 

F9.0 

29-36 

Meters 

A  priori  range 

S(l) 

F7.0 

38-45 

Arc  sec 

Azimuth  sigma 

'  ?• 

5(2) 

F7.0 

46-53 

Arc  sec 

Elevation  sigma 

1 

5(3) 

F7.0 

54-61 

Meters 

Range  sigma 

'i 

\ 

U1 

F12.0 

62-73 

Linear  error  if  IS/Q 

V 

i' 

1 

SU 

F7.0 

74-80 

Origin  error 

Vi. 

& 

£' 

2* 

XD 

F12.0 

25-36 

Degrees 

Station  latitude 

s- 

% 

£ 

YD 

F12.0 

37-48 

Degrees 

Station  longitude 

| 

ft. 

HT 

F12.0 

49-60 

Meters 

Station  height 

i, 

f 

3* 

AX 

F12.0 

1-12 

Meters 

Earth  axis 

f. 

£ 

ECC 

F12.0 

13-24 

Eccentricity  squared 

1 

XP 

F  12.0 

25-36 

Degrees 

Station  latitude 

I 

YP 

F12.0 

37-48 

Degrees 

Station  longitude 

1 

HP 

F12.0 

49-60 

Meters 

Station  height 

*  Station  position  for  the  two  baseline  stations  involved. 
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Card  No.  No  mo  Hold  Columns  Units 

If  IS  f  0 

4  COE(l)  F11.0  1-11 

*  •  • 

COE (7)  FI  1.0  67-77 


Description 


Linear  coefficients  between 
station  pair 


Repeat  the  above  cards  for  each  baseline  desired,  then  follow  the  last  baseline  set  with  a 
blank  card. 


BLANK 


Specification  Cards 


Card  No. 
1 


Name 

Field 

Columns 

Units 

Description 

NPASS 

F13.0 

1-13 

Number  of  passes 

ECC 

F13.0 

14-26 

Eccentricity  squared 

AX 

F13.0 

27-39 

Meters 

Semi-majoi  axis 

GR 

F13.0 

40-52 

Meters3 /Secs 

Gravitation  constant 

ROT 

F13.0 

53-65 

Radians/Sec 

Rotation  rate  of  earth 

IDNT 

17 

1-  7 

Station  ID 

Gl 

F7.0 

8-14 

Degrees 

Latitude 

G2 

F7.0 

15-21 

Minutes 

G3 

F7.0 

22-28 

Seconds 

G4 

F7.0 

29-35 

Degrees 

Longitude  (West) 

G5 

F7.0 

36-42 

Minutes 

G6 

F7.0 

43-49 

Seconds 

G7 

F7.0 

50-56 

Meters 

Height  (above  spheroi 

CM  (1) 

F7.0 

57-63 

Seconds 

Sigma  latitude 

CM  (2) 

F7.0 

64-71 

Seconds 

Sigma  longitude 

CM  (3) 

F7.0 

72-78 

Meters 

Sigma  height 

Card  No.  Name  Field  Columns 


Description 


2B 

23  is  the 

same  as 

2A  for  the  next  station. 

These  are  station  cards.  Follow 

the  last  station  card  with  a 

blank. 

2C 

Blank 

■ 

2D 

Ox 

F13.0 

1-13 

Meters 

Sigma  of  center  of  mass 

Oy 

F13.0 

14-26 

Meters 

Or 

FI3.0 

27-39 

Meters 

3A 

ST  11(1) 

Fil.O 

1-11 

Meters 

Input  sigma  of  plate  x 

2 

12-22 

Unitless 

Correlation  coefficient 

3 

23-33 

Meters 

Sigma  of  plate  y 

4 

34-44 

Unitless 

Correlation  coefficient 

5 

45-55 

Seconds 

Sigma  of  time 

6 

56-66 

Meters 

Sigma  of  range 

7 

67-77 

Unitless 

Correlation  coefficient 

3B&3C  38  and  3C  are  the  second  and  third  alternate  schedules  of  error  inputs  for 

plate  or  pass  data.  3A  Is  the  standard  schedule  of  error  input.  Cards  3  and 
4  are  for  optical  observations. 


4A 

ST21  (1)  F13.0 

1-13 

Radians 

j 

1 

2 

14-26 

Radians  j 

>  Orientation  sigma  a,  X 

3 

27-39 

Radians 

1 

4 

40-52 

Meters 

Focal  length  sigma 

5 

53-65 

Seconds 

Interstation  timing  sigma 

4B&4C 

4B  and  4C  are  the 

alternate 

schedules  of  sigmas  for  station  measurements. 

ST31  (1 )  F13.0 

1-13 

Meters 

2 

14-26 

Meters 

3 

27-39 

Meters 

4 

40-52 

Radians 

5 

53-65 

Seconds 

6 

66-78 

Radians 
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Error  coefficients  of 
electronic  error  model 


Card  No.  Nome  Field  Columns  Un i ts  Description 

5B&5C  53  and  5C  are  alternate  schedules  for  range  observations. 


ST41  (1)  F13.0  1-13  Meters  I 

2  1 4-26  Meters  1 

3  27-39  Meters  J 

4  40-52  Meters/Sec 

5  53-65  Meters/Scc 

6  66-78  Meters/Sec 


Sigmas  of  orbital  initial 
conditions  (position) 


Sigmas  of  orbital  initial 
conditions  (velocity) 


6B&6C  6B  and  6C  are  alternate  schedules  of  orbit  sigmas. 


Orbit  Input 

Each  orbit  will  be  set  up  as  shown  below. 


SI 

7.0 

Hours 

S2  ' 

7.0 

Minutes 

S3 

12.0 

Seconds 

Cl 

-  7.0 

Hours 

C2 

7.0 

Minutes 

C3 

12.0 

Seconds 

PI 

7.0 

hours 

P2 

7.0 

Minutes 

P3 

12.0 

Seconds 

Time  of  initial  conditions 


Time  of  desired  epoch 

Hour  angle  of  Greenwich  forzero 
hour  of  day  of  initial  conditions 
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Card  No. 

Name 

Field 

Columns 

Units 

Description 

2 

ilN(l) 

F13.0 

1-13 

Meters 

Initial  position 

2 

14-26 

3 

27-39 

4 

40-52 

Meters/Sec 

Initial  velocity 

5 

53-65 

6 

66-78 

3 

IDPAS 

15 

1-  5 

Pass  identification 

IORB 

15 

6-10 

Selects  orbit  sigma  schedule 

NSTA 

15 

11-15 

Number  of  stations  in  pass 

NINT 

15 

16-20 

Number  of  intervals  observed 

4* 

FLEA 

F12.0 

1-12 

Meters 

Focal  length  if  optical 

ITYPE 

18 

13-20 

-  is  electronic/  +  is  optical 

IS 

18 

21-28 

Selects  station  sigma  schedule 

IB 

18 

29-36 

Selects  optical  plate  schedule 

NP 

18 

37-44 

Number  of  intervals  observed 

NSD 

18 

45-52 

Identification  of  interval  observed 

*  Each  station  observing  will  have  a  card  defined  by  4.  These  cards  define 
all  data  information  for  the  particular  station. 

5 

INT 

14 

1-  4 

Number  of  total  intervals  observed 

HR 

F10.4 

5-14 

Hours 

Start  time  of  observations 

MIN 

F10.4 

15-24 

Minutes 

SEC 

F10.4 

25-34 

Seconds 

XV 

F10.4 

35-44 

Hours 

Stop  time  of  observations 

XN 

FI  0.4 

45-54 

Minutes 

XM 

F10.4 

55-64 

Seconds 
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The  following  ccjkIs  are  observation  cards  and  will  follow  the  same  order  as  the 
station  data  definition  cards  4.  Let  6  represent  an  op  ica[  obsc;vation  set  and  7 
electronic.  A  blank  card  follows  each  set  of  observations. 


Card  No. 

Name 

Field 

Columns 

Units 

Description 

6A 

GEOS  formatted  data  in  right 

ascension  and  declination 

7  A 

GEOS  formatted  data  in  range 

6B 

IID 

17  ‘ 

1-  7 

Station  Identification  number 

TO) 

F13.0 

8-20 

Seconds 

Time  of  observation 

X(l) 

F13.0 

24-28 

Meters 

Plate  X  coordinate 

Y(l) 

F13.0 

29-46 

Seconds 

Plate  Y  coordinate 

7B 

IID 

17 

1-  7 

Station  identification  number 

TO) 

F17.0 

8-24 

Seconds 

Time  of  observation 

RANGE 

F17.0 

25-41 

Meters 

Range  measurement 

B.  Output 

Program  PREP  outputs  the  primary  input  control  parameters  which  define 
the  error  parameters  and  survey.  This  provides  the  user  with  later  references  to 
the  network  characteristics,  participating  stations  and  a  priori  errors  of  the 
particular  computer  run.  Other  output  is  adjusted  parameters,  residuals, 
standard  deviations  of  recovered  parameters,  correction  to  all  adjustable  parameters 
and  total  corrections  to  survey  parameters. 

For  each  orbit  in  the  solution  the  following  output  is  written. 

1)  Residuals  for  each  observation  for  each  station. 

2)  Sigma  (standard  aeviation)  of  recovered  timing  and  orientation  biases 

for  optical  measurements. 
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FIGURE  2.  DECK  SETUP 
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6.0  GLOSSARY  OF  TERMS 


Program  Constants 

Program  Name 

Quantity 

Description 

RAD 

.01745329252 

Radian  to  degrees  conversion 

GR 

Input(optional) 

Gravitation  constant 

ROT 

Input(optional) 

Earth  rotation  rote  (mefers^/seconcP) 

AXIS 

In  put  (optional) 

Earth  semi-major  axis 

ECC 

\ 

Input(optional) 

Earth  eccentricity  squared 

CNM  ( 

SNM  j 

Smithsonian  standard  earth  model  for  1 966  through  M,  N  =  4,4. 

Symbol  Correlation 

Program  Name 

Math  Symbol 

Description 

ROT 

* 

Earth  rotation  rate 

GR 

M. 

/ 

Earth  gravitation  constant 

STD(l) 

<Jx 

Standard  error  in  plate  x 

STD  (2) 

Px 

Correlation  coefficient  in  plate  x 

STD  (3) 

ay 

Standard  error  in  plate  y 

STD  (4) 

Py 

Correlation  coefficient  in  plate  y 

STD  (5) 

°T 

Standard  error  in  timing 

STD  {6) 

Pr 

Standard  error  in  range 

STD  (7) 

Pr 

Correlation  coefficient  of  range 

STER(l) 

„  i 

a  u 

) 

SLr.R(2) 

j 

Error  in  orientation 

STER<3) 

Px  1 

) 

STER(4) 

Pc 

Error  in  focal  length 

STEP.  (5) 

Pt 

Interstation  timing  bias 

STEL(l) 

Zero  set  error  constraint 

STEL(2) 

% 

Interstation  timing  bias  error  constraint 

STEL(3) 

% 

Frequency  bias  constraint  (satellite oscillator) 

ST  EL  (4) 

% 

Frequency  bias  constraint  (ground oscillator) 

ST  E  1.(5) 

O’. 

4 

Frequency  drift  constraint 

ST  EL  (6) 

^6 

-33- 

Residual  refraction  error  constraint 

Prop  ran  Name 

Math  Symbol 

Doscriotion 

ORB(l) 

\ 

ORB  (2) 

u  1 

‘  7  X 

ORB(3) 

a  f 

0RB(4) 

°i  ( 

Standard  error  of  position  and  velocity 

ORB(5) 

<T  -  1 

y  1 

ORB  (6) 

^  l 

R1 

R; 

Rotation  matrix  (inertial  to  earth  fixed) 

XE 

xt  ) 

YE 

Orbit  position  for  the  j  th  observation 

YJ  ( 

(earth  fixed) 

ZE 

z3  ) 

Earth  fixed  coordinates  of  the  I  th 
station 


XP 

*r ) 

>  Computed  plate  coordinates 

■? 

j 

i 

| 

YP 

y°;°  j 

I 

•s 

XJ 

*°>  } 
1 

1 

>  Measured  plate  coordinates 

-a 

Ji 

YJ 

y°  ' 

> 

A 

TJ 

Time  of  j  th  measurement 

T 

$ 

RJ 

«°i° 

Computed  range 

\ 

XJ 

Measured  range 

3 

1 

DR 

Refraction  correction  to  range 

% 

'i 

WDD 

K 

A  priori  weight  matrix 

-n 

I 

WWW* 


Progiom  Name  math  Symboi  V  ,‘Scription 


s 

A 

6 

Station  correction  vector 

EPO 

5w 

Orbit  correction  vector 

EPP 

«» 

Correction  vector  to  error  parameters 

NSTA 

'Ac 

Number  of  stations  in  kth  pass 

NPASS 

K 

Number  of  passes  in  network 

1C 

J 

Number  of  observations 

XDI 

X  ' 
1 

1 

YDI 

Y 

►  Earth  fixed  velocity 

ZDI 

.  i 
Z  t 

) 

VK 

o 

Coefficients  of  matrizant  polynomials 
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GENERAL  DATA  FLOW 


I 


: ; 


Pro  Processor 

A.  General  Description 

The  Pre  Processor  (PREPRO)  program  was  designed  as  a  supplement  to  SAGA.  !t 
makes  corrections  to  data  for  polar  motion,  parallactic  refraction,  time  corrections 
(UTC-UT1  or  SA0-UT1),  and  phase  angle.  At  option  the  data  is  fit  to  a  series  of 
polynomials  (2nd  thru  6th  order)  for  random  measurement  residuals  evaluation. 

In  addition  to  optical  corrections,  PREPRO  also  converts  GEOCEIVER  observations 
(Doppler  counts)  to  range  differences  corrected  for  ionospheric  refraction.  The  residual 
option  also  exist  for  range  measurements. 

Care  must  be  taken  in  the  application  of  the  admissible  preprocessor  corrections 
so  that  corrections  that  have  already  been  made  are  not  duplicated.  Thus  the  fact  thot 
PC-1000  data  processed  to  date  by  ACIC  have  not  been  corrected  for  polar  motion, 
parallactic  refraction  or  for  UTC  to  UT1  does  not  preclude  that  at  some  time  in  the  future 
ACIC  policy  may  change  in  this  regard.  When  and  if  it  does,  the  corrections  should  no 
longer  be  applied  in  the  Pre  Processor.  Thus  one  should  remain  up  to  date  on  the  policies 
of  the  various  agencies  that  provide  data.  In  the  table  below  it  is  indicated  which 
corrections  should  be  applied  in  the  Pre  Processor  to  the  various  types  of  data  as  of 
January  1969.  Though  not  listed  in  the  table,-  corrections  for  phase  of  optically  observed 
spherical  passive  satellites  can  also  be  generated  on  option  by  the  Pre  Processor. 
Corrections,  when  needed,  for  tropospheric  ranging  refraction  are  applied  in  the  main 
program  rather  than  in  ihe  Pre  Processor. 
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T ho  user  selects  from  Table  I  the  corrections  desired.  Non-zero  values  are  inserted 
.  the  first  card  corresponding  to  the  particular  correction.  For  each  correction  desired 
the  foliow-up  inpur  cards  furnish  needed  input  paiamcters  to  accomplish  the  correction. 
These  ore  outlined  according  to  each  correction. 

Specification  Card 


Name 

Field 

Columns  Units 

Description 

ITPE 

14 

1-  4 

Type  of  measurement  (Range  or 
optic) 

'POL 

14 

5-  8 

Polar  motion 

i  REF 

14 

9-12 

Parallactic  refraction 

ICOR 

14 

13-16 

UTC  to  UT1  time 

IRES 

14 

17-20 

Residual  computation 

1  PH  S 

14 

21-24 

Phase  angle  corrects  n 

IS  AO 

14 

25-28 

ITPE  -  1 

Optical  Data 

SAO  to  UT1  correction 

H 

F4.0 

1-  4 

Degrees 

Station  •atitude 

B 

F4.0 

5-  8 

Minutes 

D 

F7.0 

9-15 

Seconds 

T 

F7.0 

16-22 

Degrees 

Station  longitude 

Z 

F7.0 

28-29 

Minutes 

V 

F7.0 

30-36 

Seconds 

HE 

F7.0 

37-43 

Meters 

Station  height 

F 

F9.0 

44-52 

Meters 

Focal  !ength 

IDD 

14 

53-56 

Station  identification 
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Card  Mo. 

K^mo 

Mold 

Columns 

Units 

Description 

f  ^ 

]t 

•j 

•  i 

F4.0 

1-  A 

Degrees 

Greenwich  hour  angle 

o 

14.0 

5-  8 

Minutes 

For  day  of  interest 

| 

k 

s 

F7.0 

9-15 

Seconus 

! 

1 

» 

SHT 

F7.0 

16-77 

Naulical  Miles 

Satellite  height 

* 

? 

v; 

COK 

F7.0 

7.3-29 

Seconds 

UTC-Ufl  correction 

I 

XAN 

F7.0 

30-36 

Arc  seconds 

x  polar  angle 

YAN 

F7.0 

37-43 

Arc  seconds 

y  polar  angle 

i 

SCAL 

F9.0 

44-52 

Scale  on  residuals 

! 

Insert  Only  If  ISAO  /  0 

!  4A 

( 

TAU 

F4.0 

1-  4 

Years 

Fraction  of  tropical  year 

1 

1 

TZ 

F4.0 

5-  8 

Epoch  of  equinox 

i 

i 

j 

DEL 

F7.0 

9-15 

Radians 

Nutation  in  longitude 

OSLO 

F7.0 

16-22 

Radians 

Obliquity  of  the  ecliptic 

i 

GTD 

F7.0 

23-29 

Seconds 

A1  time  minus  UTC  time 

i 

f 

Insert  Only 

r  IPHS  7  0 

c 

i  /n 

-  4C> 

RADI 

FI  0.0 

1-10 

k  A  L  .  .. 
iVIC  1  Cl  5 

Radius  of  satellite 

’ 

HEG. 

.  10.0 

11-20 

Meters 

Height  of  satellite 

j 

RAS 

F 1 0.0 

21-30 

Degrees 

Right  ascension  of  sun 

DES 

FI  0.0 

31-40 

Degrees 

Declination  of  sun 

5  DATA  -  GEOS  Format 


6  Blank 


i 

I 
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Card  No. 

Name 

Fie  id 

Columns  Units 

Description 

ITPE  2 

Range  Data 

Fol low  s p. 

ccification 

card  with  GEOS 

formatted  range  data. 

Follow  the  data  with  a  blank 

card  to  to 

■rminate  the 

set . 

ITPE  =  3 

GEOCEIVER  Data 

2 

cs 

F10.0 

1-10 

Mts/sec 

Light  speed 

FO 

F10.0 

11-20 

Cyc/sec 

GEOCEIVER  reference  frequency 

FS 

F10.0 

21-30 

Cyc/sec 

Transmitted  frequency 

3 

ID 

no 

1-10 

Station  identification 

X(!) 

F18.9 

11-28 

Time 

Y(l) 

F 1 8. 9 

29-47 

Doppler  count  measurement 

RR 

FI  8.9 

48-65 

Refractive  Doppler  count 

4 

Blank 

C.  Output 

The  output  is  in  punched  card  form.  The  format  is  compatible  with  SAGA  input.  The 
only  written  output  is  the  residuals  of  the  2nd  thru  6th  order  polynomials  and  conic. 

Written  output  of  the  corrected  observations  may  be  obtained  by  changing  the  file 
unit  number  to  the  appropriate  output  file  number. 

D .  Analysis 

Polar  Motion  Correction 

Let; 

T  =  time  of  observation  in  years  and  days  (e.g.,  T=  1966.385) 

t  =  time  of  observation  from  ZULU  midnight 

rjta-  right  ascension  observed 

x,y=  angles  of  polar  motion  (radians). 
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^  o.n  pL: !  s.‘ : 

(1)  S  -  apparent  Greenwich  sidereal  time  cil  lime  r  of  observation 

=  P0  +  an,  whore  9-,  =  Greenwich  hour  angle 
cc-  -  earth's  rotational  rate 

(2)  In  terms  of  given  ey,  5  .‘r,d  computed  value  of  0  <,y  evaluate  thr  expression 


r 

\ 

' 

f 

■ 

X 

cos  0  -sin  0  0 

1  0  x 

cos  0  sin  0  0 

cos  ry  cos  6 

sin  6  cos  0  0 

0  1  -y 

' 

-sin  0  cos  0  0 

sin  o  cos  ft 

V 

0  0  1 

-x  y 

0  0  1 

sin  ft 

* 

, 

L 

. 

„  > 

(3)  Compute  r/ ,  ft'  from  the  relations 

sin  <v'  =  fj  /  ~v' 

cos  r/  -  X/  n’1  ~v-‘- 

sin  ft1  =  v  (sign  of  ft'  is  same  as  sign  of  v) 

(4)  Replace  ry,  ft  by  o' ,  6' . 

Phase  Angle  Correction 

The  apparent  direction  of  a  spherical,  sun  reflocihig,  satellite  may  be  biased, 
depending  upon  the  portion  of  surface  observed.  Two  different  types  of  light  reflections 
are  encountered.  First  is  that  of  diffused  reflection,  which  it  the  result  of  rejections 
from  irregular  surfaces.  The  moon  is  a  example  of  diffused  reflection.  Second  is  that 
of  specular  reflections,  the  result  of  light  reflected  from  smooth  surfaces  such  as  a  mirror. 

In  the  case  oi  n  on-spherical  sale S I •  -  .s,  the  correction  is  treated  as  a  constant  bias 
and  y  (sect ‘on  4).  An  attempt  to  apply  direction  corrections  would  be  nearly  impossible 
and  ccrtoinl)  impractical  because  of  vuriout  satellite  configurations  and  attitudes . 
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v.:  .j,  UA11.1VI1U  .4  1110  Ks, 


... ser  v«-i, iui.  iii  required  in  determining  the  direction  of 
r..o  sun  fro...  ; he  eu.th's  center.  The  time  of  either  the  pre-  or  post-calibration  is  suf¬ 
ficient.  !',  doth  are  available  the  mean  of  the  two  may  be  taken  as  the  time  of  observation. 


On.  use.':  ;C  i 1  act  ion 


Compute  the  direction  cosines  (inertial)  of  the  sun  for  time  of  observation,  given 
right  ascension  (tv),  declination  (ft)  and  time. 


X  --  cos  n  cos  t, 
(1)  u  -  sin  rt  cos  6 
v  =  sin  ?>. 


Convert  to  earth-fixed  direction  cosines 


(2) 


*  1  fees  $. 


sin  ts- 


cos  9  -sin  9, 

£  £ 


whore 


(3)  S4  =  vOo) 


6.  - 


right  ascension  of  Greenwich  for  year  of  interest 

-  rotation  rote  of  earth  in  radians/day 

-  time  in  days  and  decimal  days  from  January  0.  * 


Let 


(-')  V.  XshH  u,  i1 


be  the  direction  vector  from  the  earth’s  center  to  the  sun.  Assume  Vi  is  parallel  to  the 
vector  from  the  satellite  !o  the  sun. 
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-i  Mr»fc. rfctajd! 


from  ct ,  ft  to  satellite  and  let 


Compute  \  jj  ,  v  ( . 

(5)  V:1J  -  |_+  Vtjk 


bo  the  j  fK  dircc'ion  vector  from  station  i  fo  the  satellite. 


The  dot  product  of  (+V.)  and  (~Vr.)  yields  the  intersection  angle  at  the  satellite 


(6)  mu  =-  cos  *  ((+V,)(-Vnui) 


The  distance  dSJ  normal  to  V.-,SJ  from  the  satellire  center  is  inversely  proportional 
to  the  percentage  of  the  illuminated  surface  viewed  by  camera  i.  Solve  for  d1}: 


(7) 


d'n  - 


d‘.3 

sin  r),  j 


where  r5  is  the  sate! life  radius  and  d) ,  has  the  same  direction  as  -V x .  The  correction 


vector  Vli  i  is 


(S)  Vl4J  -  d;.A.sj_  +  diju,  1+  dJjV.k 

The  magnitude  cf  V3l,  is: 

i 

(9)  R. ,  =  |(X“  -x;f  +  +  (Z°r2‘)3j  ' . 


Then  V„  is: 


(10)  V,,.  - 


The  vector  V..  from  station  i  to  the  i  fh  sa.  ollite  center  is: 


01)  v.,;  =  v11J+vSj>. 
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Oi-OCriVER  Range  HiiTcivnco  Computation 

Sec  Section  3.0  of  i  art  I  for  a  detailed  analysis  of  the  GEOCEIVER  range 
difference  reductions  and  the  formation  of  the  basic  error  model. 
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APPENDIX  B 
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j^iin.niy  U-Cr...  .c  |;,p lection 


,  >  .*11 

v-'ji i^Oi  ojjci v of i ons 


Qie  converted  from  right  ascension  and  declination  to  plate 
'tMX'>V'  h‘Sl  COnvert  r'3ht  asconsi°n  and  declination  to  azimuth  (A,)  and 
'■  ^  '}'  A"b',T'c  tho  midd,c  observation  (azimuth  and  elevation)  to  be  the  - 


yj;  , ,\£  principal 
\'J  /X ,  x}  oecon'e; 


a^is  ot  tne  camera.  The  external  orientation  elements 


wr.e  re 


n1  '  A.  ,  at  -  E„  and  let  X  =  0. 

For  -he  ,  rb  station  end  time  compute  the  local  hour  angle  of  the  j  th  observat 
A  -  tu  +  1  .00273791  t,-X-n< 

ru  Greenwich  hour  angle  at  t  =  0  (apparent) 

.  .  -  oosorving  time  (universal)  in  angular  measurement  from  Greenwich  (radians) 

X  ~  longii'ude  of  station  (West) 

fr  :lu,H  ascension  of  observation 
''  -  taritucle  of  station 


sm  r.  =-  sm  sin  6+  cos  *  cos  6  cos  (ts) 
cos  E  =  (1  -sin  E)® 

sir,  A  -  -cos  6  sin  (tj/cos  E 

Cv.s  r,  -  (sin  6-  sin  Vs  in  E)/(cos  <r>  cos  £) 


ion . 


c  . 


tar,  (sin  A/co;  A) 
ran  (sin  E/cos  £} 
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;  Sr’4  -  '‘Y~  * 


lo  compute  the  local  orientation  matrix. 


use  midpoint  A  , ,  E  , 


f  A 

8 

f  -cos  rv 

sin 

1  A 

A 

t  f 

1  A' 

D ' 

C‘ 

1  | 

!  1 

-sin  rv  sin 

CO 

-cos 

sin 

(0 

A 

'  t 

V  D 

E 

F  J 

’ 

J 

1  sin  rv  cos 

£.0 

COS 

cos 

a.' 

then 


sin  A 
cosA 
sin  E 


cos  E 
cos  E 


the  plate  coordinates  for  the  j  th  point  are: 

=  cuj/wj 
y3  =  cv  j /w  r 

.A  W  A 

Let  (x b , y-t, )  denote  the  coordinates  of  the  first  and  last  points  on  the  trace.  Let 

X  denote  the  angle  between  the  line  joining  these  two  points  and  the  x  axis.  Then 
compute: 

sin  x  -  -(yb-yJ/Vttb 
cos  x  =  (xb  -xu)/rmll 

1 

hi  =  [(xb  -xj3  +  (yb -yj3]' 

The  new  x,  y  coordinates  are: 
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r*l 

V 


1 


-cos  X  sin  x 


sin  x  cos  x 


r- 

x 


The  rotation  of  the  local  orientation  into  the  master  frame  is  accomplished  by: 


A 

B 

C 

A1 

B' 

C' 

= 

D 

E 

F 

-cos  X  sin  x  0 

sin  x  cos  x  0 

0  0  1 


ABC 

A  A  A 

A'  B1  C‘ 

AAA 

D  E  F 


0  1  0 

-sin  <p  0  cos  (p 

cos  <p  0  sin 


cos 


sin  X  0 


sin  X  cos  X  0 


0  1 


The  coordinates  x^y^  and  the  master  orientation  matrix  serve  as  input  to  the  main  program. 
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Interstation  Constraints 

A.  Incorporation  of  Special  Infestation  Constraints 

From  the  latest  approximations  to  the  coordinates  of  the  j  fh  specified  station  pair 
pj,  q  ,  compute  (subscripts  j  are  omitted  in  following): 


AXp, 

=  (>4)00 

-  «)00 

0) 

=  (Yp)00 

-  (Y^)00 

“  (Z'p)00 

-  (2?/° 

R°°  = 

C(*xp,)! 

3  +  (^YPJ2  +(^ZM)3] 

C(Axn)s 

+  (AY^)a  £ 

sin  Ap°  =  .jYM/r°°  ,  A°°  =  arc  sin  (sin  A°°) 
sin  E*  =  AZk/R°°7  E°m°  «  arc  sin  (sin  E* ) 

^xw/r“ 


t?  i 

* 


1  &A, 


I  cosA_ 


i  aE. 


1  (X„sinE,3 


<T*  w 

3Xp 

rw 

bu  = 

ar 

BXj, 

»q 

rp<i 

1 

SA00 

1 

sinA°N° 

l  - 

1 

< 

1  (*: 

•  c00\ 
'pqs,n  Epq) 

<>am 

w, 

SYP 

crA 

an 

rw 

'V4 

crE 

eM 

BYP 

Tpij 

1 

< 

■  =  0 

U  - 

1 

BE* 

1  0  " 

«'**,nES) 

Q 

> 

S 

3ZP 

Du  - 

azp 

<7epi 

*'P4 

i  »”  i 
Cu  = «. 


cia  " 


1  BRp,  1 

wr=^r~il^ 

w  si 


l  aR00  i 

=  5TV'W' 


When  a  prespecified  linear  constrainf  is  to  be  imposed  on  the  coordinates  of  stations 
p  qnd  q  compute: 

(3)  U°p°  =<**  (^p)00+aa(Y£)00  +Oi9(Zf:J))00+  j31(Xp)“  +  ^a(Vp00+|59(2^)00 


i  au00  i 

dll=ai>T=a“l 


i  au00  i 
axf=; 


,  _  i  au00  i 

dia-  “  -^"  =  -  “a 
cr  oYp  cr 


d  -1  W°°  -1  ft 

16 %  SyT'o  ®* 


.  .  I  3U®  1 

di,‘;  izf"?  “■ 


,  _  i  au00  i 
d*~*  azf  =  ; 


72= 


In  terms  of  the  above  sot  up: 


<■*> 


an 

alo 

aia 

-Oji 

“Oja 

-al3 

bn 

bJ3 

~bu 

“bia 

-bln 

cn 

cia 

“cli 

“C1B 

“C14 

Jn 

dla 

“du 

“dia 

-d:s 

Azimuth  Constraint 
Elevation  Constraint 
Distance  Constraint 
Linear  Constraint 


Let  Ap,, denote  the  measured  values  and  let  <7*  P<1  »  ^tp,  or*p(i  crUp^ 
be  the  corresponding  standard  deviations.  Then  set  up  the  discrepancy  vector: 


<-A°pq)/crAp<i 

C-U°pq)/( 


If  a  particular  type  of  constraint  is  not  to  be  exercised  between  points  p  and  q,  the  rows 
of  Up,  and  ep,  corresponding  to  that  constraint  should  be  set  equal  to  zero.  Thus  if  a 
distance  constraint  were  to  be  rows  of  Up,  end  cp,  would  consist  of  zero  elements.  In 
terms  of  Up,  and  ep,  compute: 


(SX6) 


^pq  ^pq'  cpq  ^pq  €pq 


Partition  Sp,and  cp,  as  follows: 


S  S 

jpP  -'p, 


(a,a)  (s,a) 


. T  . 

S  S 

Jqq 

(a,a)  (a,®) 


f  cpq 
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(6)  Merge  Sp?  and  cp4  info  the  proper  position  of  the  network  normal  equation,-  S  j 

$ 

and  c,  respectively.  1 

I 

(7)  Continue  as  above  until  all  intorstation  constraints  have  thus  been  processed.  | 


i 
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APPENDIX  D 


Preliminary  Geoceiver  Computations 

The  raw  Geoceiver  cycle  counts  AN  j  are  first  converted  to  range  differences 
between  the  satellite  over  the  time  interval  tJ-xto  tj.  If 

f0  =  transmitted  frequency 

fg  «  Geoceiver  reference  frequency 

tj  -  time  of  j  th  cycle  count 

Sj  =  slant  range  between  Geoceiver  and  satellite 

c  =  velocity  of  light 

then 

=  *  j  ”  ~  tfo  “fo)fr  j-i) 

where 

\>  "  «A> 

and  Dj  is  the  cycle  count  corrected  for  ionospheric  refraction.  Dj  is  given  by: 

Dj  =  ANj-  K  AAN 3 

where 

AANj  =  refractive  cycle  count 
K  =9  for  162-324mc  reception 

K  =9  1/6  for  15Q-400mc  reception. 

The  range  differences  are  converted  to  relative  ranges  by: 

rj  ~  Asj  +  Asa  +  . . .  +  Asj 

where 

r3  =  change  in  range  from  time  t  =  t0  to  time  t  =  t3 . 

The  quantities  t3  ,r3  constitute  the  input  to  the  main  program.  When  any  particular 
AN3  is  equal  to  zero,  this  means  that  phase  lock  was  lost  during  the  jth  counting  interval 
and  it  becomes  necessary  in  the  main  program  to  reinitialize  zero  set  at  time  t }  . 
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